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SUMMARY 


A semi-elliptic formulation, termed the interacting parabolized Navier- 
Stokes (IPNS) formulation, is developed for the analysis of a class of 
subsonic viscous flows for which streamwise diffusion is negligible but 
which are significantly influenced by upstream interactions. The IPNS 
equations are obtained from the Navier-Stokes equations by dropping the 
streamwise viscous-diffusion terms but retaining upstream influence via the 
streamwise pressure-gradient. A two-step alternating-direction-explicit 
numerical scheme is developed to solve these equations. The quasi- 
linearization and discretization of the equations are carefully examined so 
that no artificial viscosity is added externally to the scheme. Also, 
solutions to compressible as well as nearly incompressible flows are 
obtained without any modification either in the analysis or in the solution 

procedure. 

The procedure is applied to constricted channels and cascade passages 
formed by airfoils of various shapes. These geometries are represented 
using numerically generated general curvilinear boundary-oriented 
coordinates forming an H-grid. Stagnation pressure, stagnation temperature 
and streamline slope are prescribed at inflow, while static pressure is 
prescribed at the outflow boundary. Results are obtained for various values 
of Reynolds number, thickness ratio and Mach number. The regular behavior 
of the solutions demonstrates that the technique is viable for flows with 
strong interactions, arising due to either boundary-layer separation or the 
presence of sharp leading/trailing edges. Mesh refinement studies are 
conducted to verify the accuracy of the results obtained. 


’RECEDING PAGE BLANK NOT FILtotL 





A new hybrid C-H grid, more appropriate for cascades of airfoils with 
rounded leading edges, is also developed. Appropriate decomposition of the 
physical domain leads to a multi-block computational domain bounded only by 
the physical-problem boundaries. This permits development of a composite 
solution procedure which, unlike most found in literature, is not a patching 
procedure. Satisfactory results are obtained for flows through cascades of 
Joukowski airfoils. The implementation of the IPNS formulation on the C-H 
grid exposes two small portions of the grid interfaces and these require 
special treatment. However, with a hybrid grid, the use of complete Navier- 
Stokes equations is recommended, so as also to avoid inconsistencies in the 
parabolization approximation due to changing orientation of the coordinates 
at a given location. 
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NOMENCLATURE 
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E , F 
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h 

i.j 
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J 


L 
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p 

Pr 


P 


0 


q 
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R 

Re 
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viii 



t 

T 


u,v 


U.V 


U 
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Temperature 

Total temperature 

Cartesian components of velocity 

Contravariant velocity components 

Mass~averaged velocity at the inflow boundary 

Cartesian coordinates 

Airfoil plane in the Joukowski transformation 


Greek Symbols 
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A 
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At 
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e abs* e rel 
C 

0 

P 
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Description 

Two-dimensional Laplacian operator 
Ratio of specific heats 
Forward difference operator 

Incremental mesh spacings in the Cartesian coordinates 
Time step 

Incremental mesh spacings in the transformed coordinate 
system 

Absolute and relative errors 

Circle plane in the Joukowski transformation 

Flow direction at the inflow boundary 

Coefficient of viscosity 

Transformed coordinates 

Metric derivatives 

Density 

Stress component, transformed time 


IX 


Wall shear parameter 

Parameter for splitting the pressure gradient term 


Subscripts Description 

b Body, backward difference 

f Forward difference 

re f Reference quantity 

v Viscous term 

w Wall 

x,y Partial derivatives wth respect to Cartesian 

coordinates 

^ ,n Partial derivatives with respect to the transformed 

coordinates 


Supers cripts Description 

i Streamwise mesh index 

j Normal mesh index 

n Temporal index 

* Dimensional quantity 
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CHAPTER 1 


INTRODUCTION 

The flow through compressors and turbines of gas-turbine engines is 
fairly complex. The complexities arise due to unsteadiness, separation, 
periodic transition from laminar to turbulent flows and complex 
geometries. A clear understanding of these flow phenomena is needed in 
order to improve the performance of these components of the engine. It 
is well known that the complete Navier-Stokes (NS) equations accurately 
describe the important physical aspects of fluid flow occurring in these 
components. However, in spite of all the advances made to date in 
numerical algorithms and computer firmware, numerical solution of the 
complete NS equations can still require large amounts of computer 
resources in terms of time and storage. Hence, an approximate form of 
the NS equations which accurately depicts the physics is preferred. The 
simplest of the approximate forms of the NS equations is provided by the 
boundary-layer equations. The classical boundary layer (CBL) equations 
with specified pressure gradient are parabolic in nature. Therefore, a 
spatial-marching procedure can be employed to numerically solve these 
equations very efficiently. However, this formulation does not contain 
any mechanism for transmitting downstream disturbances upstream and, 
hence, cannot be employed for problems where there is a strong pressure 
interaction or when the flow is separated. Goldstein [1] showed that 
the solutions to the classical boundary-layer equations exhibit a square 
root singularity in the wall shear at the separation point. This 
singularity leads to the failure of the weak interaction method wherein 


the outer inviscid flow and the inner boundary-layer region are analyzed 
sequentially, with the interaction between the two regions being 
modelled through the pressure gradient term. These limitations were 
overcome by the development of interact ing-boundary layer [2] and triple 
deck [3] theories. In interacting boundary-layer (IBL) theory the 
pressure gradient is treated as unknown. In subsonic flows, the 
pressure gradient is related to the derivative of the displacement 
thickness through Cauchy’s integral. Detailed discussion on interacting 
boundary— layer theory has been given by Veldman C2]. An interacting 
boundary-layer model has been used by Rothmayer [J»] for analyzing high 
Reynolds number flows with large regions of separation. However, the 
interacting boundary-layer model also has its drawbacks. For complex 
flows, relating the pressure gradient to the displacement thickness is 
not sufficient. Also, for flow past bodies with large curvature, the 
normal pressure gradient is no longer negligible and should be included. 

To account for these effects, Briley [ 5 ] and Ghia et al. [6] 
developed a non-iterative parabolic procedure for calculating flow 
through curved ducts. Their procedure employed parabolized Navier- 
Stokes ( PNS ) equations obtained by neglecting the viscous diffusion 
terms in the streamwise direction, with the streamwise pressure gradient 
term being represented by a backward difference. Hence, this procedure 

is applicable for flows with little upstream influence and no streamwise 
separation. 

The thin-layer Navier-Stokes (TLNS) equations of Steger [7] 
include the upstream influence. These equations are obtained from the 
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unsteady NS equations by dropping the strearawise diffusion terms. The 
procedure employed to solve these equations is a -time-marching’ 
technique and has proved to be costly in terms of computer time, in 
order to obtain steady-state solutions of flows around isolated 
airfoils. Steger, Pulliam and Chima [8] have employed the two- 
dimensional TLNS equations and a C-type of grid for solving viscous 
flows through cascades. They experienced difficulties in obtaining 
steady-state solutions when the pressure is not prescribed at the 
upstream boundary. Buggeln, Briley and McDonald [9] have computed 
laminar and turbulent flows through ducts using the Navier-Stokes 
equations. Chima and Johnson [10] employed an explicit multiple-grid 
algorithm to solve the NS equations in order to improve convergence. 
Shamroth, McDonald and Briley [11] and Hah [12] have computed cascade 
flows using the complete NS equations. Rhie [13] bas employed 
parti ally-parabolic MS equations to analyze three-dimensional viscous 
flows through curved ducts of arbitrary cross-section. Recently, 

[lit], Davis et al. [15] and Rhie [16] have developed methods for 
predicting cascade flows using NS equations. References [U] and [16] 
have also employed a multigrid algorithm to enhance convergence. Most 
of the works mentioned above have incorporated second- and fourth-order 
dissipation terms, in order to suppress oscillations in the flow field. 
The difference in computational effort involved in obtaining the 
solution to TLNS and complete NS equations is not significant. The 
numerical solution of both the TLNS and the complete NS equations 
require large amounts of computer resources. 
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In the present study, a single system of equations which can 
include the upstream influence is obtained from the full NS equations. 
It is termed the interacting parabolized Navier-Stokes (IPNS) 
formulation and belongs to the class of semi-elliptic models, one form 
of which was developed earlier by U. Ghia et al. [17]. Only steady 
flows are discussed here and, hence, the time-derivative term in the NS 
equations is dropped. It should be mentioned, however, that the 
analysis can be extended readily to unsteady flows by the inclusion of 
this term. The serai-elliptic form of the equations is obtained by 
dropping the viscous diffusion terms in the streamwise direction. This 
approximation is supported by the fact that the streamwise diffusion is 
negligible compared to the normal diffusion, for the flows under 
consideration. Clearly, the approximation is appropriate if the 
coordinate system employed is a body-oriented, near-orthogonal system. 
The semi-elliptic or IPNS formulation is tested via application to 2-D 
flows through channels with varying cross section in the streamwise 
direction and flows through cascades of airfoil 3 of various shapes. 

These configurations are chosen as they are akin to the geometries of a 
turbomachinery compressor or turbine. 

In all of the works mentioned above, either an H- or a C-type of 
grid is employed. In order to analyze flow around airfoils with rounded 
leading edges, it is often desired to employ a combination of these 
types of grids. Near the leading edge, the channel or the H-type of 
grid becomes excessively skewed and non- orthogonal and a C-grid is more 
suitable in this region. But, in the latter, the grid density decreases 
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rapidly with distance away from the leading edge. In this region, an H- 
grid can be employed. Norton, Thompkins and Haimes [18] have employed a 
mixed sheared and O-type grid for computing flows through turbine 
cascades. Rai [19] has employed a patched and overlaid grid system in 
order to compute flow through a rotor-stator combination of a 
turbomachine. Bush [20] developed a zonal methodology and a time- 
dependent procedure to obtain solution of the NS equations for flow 
through an external compression inlet. When the zonal or overlaid grid 
systems are employed to solve the governing equations of motion, it is 
important to transfer information from one grid system to the other 
appropriately. Hence, in the present study, the hybrid C-H grid 
generation procedure developed by U. Ghia, K. Ghia and Ramamurti [21] 
for turbomachinery cascades is employed. When this hybrid C H grid is 
employed to solve the complete MS equations in a composite manner, the 
explicit transfer of information across the zonal boundaries is not 

required. 

Details of the derivation of the governing equations are given in 
Chapter 2. Also, the appropriate boundary conditions to be specified 
for solving the governing equations, for both channel and cascade 
configurations, are discussed. In Chapter 3, the numerical procedure 
employed is discussed. The appropriate form of the pressure gradient 
term and the metric terms associated with it and the implementation of 
the boundary conditions and modeling for reversed flow, are also 
included in that chapter. Results for flows through constricted 
channels and cascades of airfoils of different shapes, obtained 
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employing the channel or H~type of grid, are discussed in Chapter 4. a 
composite procedure for generating a hybrid C-H grid for cascades with 
rounded leading edges is given in Chapter 5. In Chapter 6, the 
implementation of the solution procedure for flow through a cascade of 
Joukowski airfolis using a hybrid C-H grid is discussed. Some results 
obtained are presented in this chapter. Details of the implicit 
solution of a system of equations subjected to a periodicity boundary 
condition arising in cascade flows, the discretized representation of 
the metric coefficients and the treatment of the five-sided cell 
occurring in the hybrid C-H grid are included in the appendices. 
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CHAPTER 2 


FORMULATION OF THE PROBLEM 

2.1 Basic Equations 

The governing equations for the mathematical model of fluid flow 
can be derived from the Navier-Stokes equations. The nondimensional , 
conservation form of the equations for two-dimensional laminar flow of a 
compressible fluid can be written in Cartesian coordinates as follows: 
Continuity 

If * !i <"“> * h <pv) ■ 0 (2 -' a) 

x-Momentum 


ft <pu) * h (pui * p> * h <puv) ‘ h <T xx> * h ‘V 


(2.1 b) 


y-Momentum 
3 


at (pv) * fl (puv) * fy <pv!,p) ■ h (l xy’ * ly (T yy’ 


(2.1 c) 


E nergy 

ft (pe t ) * li |(p v p)u > * h 1 <p v p)v) 


% 3 

— (ut +VT -q ) + - r — (ut + vt - q ) 
3x v xx xy 4 x 3y xy yy y 


(2.1 d) 


where p is the density, u and v are the Cartesian components of velocity 
and e is the specific total energy given in terms of specific internal 


energy e by 


u 2 + v 2 


e • e + „ 

t 2 


(2.2 a) 


The stress components and the heat flux terms can be written as 
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T 


{ ( A+2y )u + 

Xv } 

XX 

Re 

X 

y 

T 

1 

m ■ — - 

{ ( A+2y ) v + 

Xu } 

yy 

Re 

y 

X 

T 

„ I_ 

{y(u + v )} 


xy 

Re 

y x 


n a 


-U T 



n r* rp 

q x “ Re Pr (Y-1 )M 2 x 
00 


and 


q y " Re Pr (Y-1 )M 2 T y 
00 


(2.2 b-f) 


2 

According to Stokes' hypothesis, X is taken as (- — y). 

The equation of state is given by 

p - (Y-1 ) pe (2. 

The constitutive equation for viscosity is given by Sutherland's 
viscosity law 


( 1 + T ) „ 3/2 

U = — T 

(T+T) 


( 2 .*) 


where 


no°K 

T 

ref 


The Reynolds number and the Mach 

the inlet boundary and are given 

Re - (p U L)/ y _ 
ref avg H ref 


number are based on the conditions at 
as 

(2.5 a) 


and 


M 

00 


u 

avg 


/ (YR T ) 1/2 
ref 


(2.5 b) 
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l x) 



where U is the mass-averaged inflow velocity at the inlet given by 
avg 


f h 

U . d P v *5 , (2.6) 

avg rh 

o J P ds 


with h as the cascade blade spacing or the channel height, V the 
velocity normal to the inlet boundary and s the distance measured along 
the inlet boundary. 


The reference length L is the chord length of the airfoils for 
cascade flows and the channel height for channel flows. 

Equation (2.1) has been Obtained by the following nondimension- 

alizat ion: 


x 




t * 


(L / U ) 
avg 


u 


* 

u 


avg 


v 


# 

V 

U 

avg 



P - u‘ 
K ref avg 


avg 


and 


* 



All the dimensional quantities are denoted with a superscript asterisk. 
Equation (2.1) can be written in a vector form as 


w + 3E ♦ 

3t 3x 3y 




( 2 . 8 ) 


where 


Q 


C P. 


pu, pv, pe t 
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9 


E - [ pu, pu 2 +p, puv , (pe t +p)u , 

F - [ pv, puv, pv 2 +p, (pe t +p)v ] T 

E ■ [ 0, t , T , (ut +vt —a ) 1^" 
v xx’ xy’ xx xy q x' J 


and 


F V ' [ V v (UT x/ v ' t yy'V ]T • 


(2.9 a-e) 


2.2 Coordinate Transformation 

The success of a numerical solution procedure for the governing 
equations of motion depends heavily on the proper choice of coordinates. 
One of the first requirements placed on a coordinate system is that the 
coordinates be aligned with the problem boundaries. The use of 
boundary-fitted coordinates reduces the complexities otherwise 
encountered in the treatment of boundaries of arbitrary shape. Hence, 
the Navier-Stokes equations in the physical (x,y) coordinates are 
transformed to a system of computational (5*n) coordinates through the 
following general transformation: 

5 - 5(x,y) , 

n - n(x,y) 

and 1 “ 1 ‘ (2.10) 
According to Viviand [22], the transformed governing equations in 
the (5,n,t) coordinates can be written in the strong-conservation-law 
(SCL) form as follows: 
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( 2.11 ) 




3 t 


3 S _ 3 n n _ 

“ H ( T E v + J- F v } + 3^ ( J~ E v + J 1 F V ) 

where J Is the Jacobian of the transformation and is defined as 


det 


9U,n) 


3(x,y) 


X 5 y n ' h \ ' 5x ^ ’ 5y " x 


( 2 . 12 ) 


The metrics E , £ » n and n„ are determined after the mapping, 

a y x y 

given by Eq. (2.10), has been defined. The metrics are related to the 
derivatives x^, y^, etc., by the following relations. 

E x ■ J % ' S ■ - J % • 


h • "y ' J x 5 


(2.13) 


It is convenient to write Eq. (2.11) in the following form. 


5Q + + 3F = 3_ 3_ 

3t 9 £ 9n 3£ 9n U v' 


(2.1 U; 


where 


°-3 


5 

E » E + -fL F 
J J 


n n 

F ■ -r~ E 


_ 5 

E ■ E +-/F 
V J V J V 


and 


n n 

F ■ -r— E + -ft F . 
V J V J V 


(2.15 a-e) 
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An alternative compact representation of the inviscid flux vectors 
E and F is also possible and is obtained by using Eqs. (2.9 b,c) in Eqs. 
(2.15 b,c) to yield 



pU 


pV 


pull + 5 X p 

and F - — j- 

puV + n x p 


pvU + p 


pvv + n y p 


(pe + p)U 


(pe fc + p)V 


where U and V are the contravariant velocities along the 5 and n 

coordinates, respectively, and are related to the Cartesian components u 

and v by 

U - 5 u + E, v 
s x y 

and 

V - n u + n v (2.17) 

x y 

Equations (2.14) are the complete Navier-Stokes equations in SCL 
form in the general boundary-oriented (£,n) coordinates. As discussed 
in Chapter 1, certain approximations will be introduced in these 
equations so as to optimize the efficiency of their numerical solution 
and their ability to represent important physical flow phenomena 
accurately. 

2.3 Derivation of the Semi-Elliptic Form of the Governing Equations 

The time-derivative term is dropped from Eq. (2.14) because only 
steady flows are considered in the present study. The parabolized 
Navier-Stokes equations are obtained by neglecting all the streamwise 
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diffusion terms. This involves dropping the second-order derivatives 
( 3 2 /9S 2 ) and the cross derivatives ( 3 2 / 3 £ 3 n ) in the viscous 
terms. This approximation is supported by the fact that the streamwise 
diffusion is negligible compared to the normal diffusion in most of the 
regions of the flows under consideration. This approximation is 
appropriate only if the (£,n) coordinate system is a body-oriented, 
near-orthogonal coordinate system, that is, the 5 coordinate is nearly 
aligned with the streamwise direction and the n coordinate is nearly 
orthogonal to it. The reduced set of equations can be written as 
follows. 


9E + 3F 
9£ 9n 



( 2 . 18 ) 


where E, F and F y are as given in Eq. (2.15). 

It should be emphasized that the above set of equations is 
'parabolized' and not parabolic. The mathematical character of the 
system of equations (2.18) depends on the manner in which the streamwise 
pressure gradient term p^ is treated. If p^ is prescribed, as in the 

case of classical boundary-layer theory, the system is parabolic. In 
this case, a marching method can be employed to obtain the solution for 
this system. This method of solution is very efficient, but it does not 
have any mechanism for including upstream influence and is, therefore, 
not suitable for flows with separation and sudden streamwise changes in 
boundary conditions. When the p^ term is treated as unknown and forward 
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differenced, the system of equations is no longer parabolic but has an 
elliptic character. 

When the parabolized Navier-Stokes equations are solved as an 
inital-value problem, as in the case of single— sweep marching solutions, 
the ill-posedness of the equations leads to ’departure solutions* 
similar to the eigensolutions of the viscous sublayer equations proposed 
by Lighthill [23]. Vigneron, Rakich and Tannehill [24] have described a 
method for suppressing the departure solutions in their study of 
supersonic flow over delta-wings. They introduced a parameter, a j, to 
split the pressure gradient term p^ into *uj P x 1 , which was backward 

differenced and treated implicitly, and (1-ui)p x * The latter term, even 


when represented using a backward difference, led to instabilities and, 
hence, was dropped entirely. This is appropriate if the flow is 
predominantly supersonic, as in the case these authors considered, but 
not in general. These authors performed a characteristics analysis for 
the inviscid as well as the viscous limits of the equations. From the 
viscous analysis, they found that the equations are well posed for space 
marching when 


y m 2 

X 


» s u(y-if TP- • f <V lf r(M *> 4 ' 


and a) - 1 if f(M ) > 1 , 


where M x 


u_ 

a 


(2.19 a) 
(2.19 b) 
(2.19 c) 


In the present study, the pressure-gradient term is split in the 
manner described above, but is discretized so as to include upstream 
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influence in flows with strong interaction by forward differencing the 


' (1-u) 3p/3£ ' term, i.e., by using the representation 


3p 




(1-u) 


3p 

K 


f 


( 2 . 20 ) 


Here, the subscripts b and f denote backward and forward differences, 
respectively. 


2.4 Boundary Conditions 

The appropriate boundary conditions to be specified for solving the 
governing equations, described in the prevoius section, for flows 
through cascade and channel configurations are described in this 
section. These conditions are specified along the inflow and outflow 
boundaries, and lower and upper wall boundaries for channels and airfoil 
portions of the cascade passages. These boundaries are shown in Fig. 1. 
Also, for cascade flows, the periodicity of the flow variables along the 
wake boundaries is considered. A set of wall-wall boundary conditions 
is one that is imposed at the ends of a ( = constant grid-line which 
starts at a body surface and ends at the facing body surface. A set of 
wake-wake boundary conditions is one which is imposed at the ends of a 
5 « constant line which starts at a point in the wake region and ends at 
a point in the facing wake region. In the case of staggered cascades, a 
set of wall-wake boundary conditions may be needed in regions where a 
£ - constant line starts from a point along a wall and ends at a point 
in the facing wake region. 
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2.4.1 Inflow and Outflow Boundary Conditions 


For the problems considered in the present study, the flow near the 
inflow and outflow boundaries behaves in an almost inviscid manner, 
von Mises [25] has carried out a characteristics analysis for inviscid 
systems and found that for subsonic flows, all the characteristics are 
real, with two of them being positive and one negative. Using the 
counting principle of Courant and Hilbert [26] that one boundary 
condition is to be specified per entering characteristic, this requires 
that two conditions are to be specified at the inflow boundary, and one 
at the outflow boundary. 

In Ref. [27], McDonald and Briley have described a specific set of 
boundary conditions. They considered a typical duct flow proceeding 
from a large reservoir and exhausting into a plenum. The reservoir 
conditions and the plenum static pressure were known. This duct flow 
model leads to prescribing the reservoir total conditions and the plenum 
static pressure. The specified stagnation temperature and pressure 
constitute the two required inflow boundary conditions and the specified 
static pressure constitutes the one outflow boundary condition. 

In the present study, at the inflow boundary, the total pressure 
prescribed for cascade flows is that corresponding to a uniform velocity 
prof ile, ' while the stagnation temperature is taken to be constant. For 
channel flows, the conditions prescribed at the inflow boundary are the 
velocity and static temperature profiles corresponding to a fully 
developed flow in the channel. 
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In the problems considered, the outflow boundary is situated far 
downstream of the cascade of airfoils or the channel constriction, so 
that uniform static pressure is an appropriate condition. 

As the static pressure at the inlet is not specified, the mass flow 
in the configuration is not set a priori and pressure waves can escape 
upstream, avoiding the problem of reflecting waves discussed by Rudy and 
Strikwerda [28]. To facilitate a marching procedure, the conditions at 
the inflow boundary are obtained by assuming the velocity-profile shape, 
guessing a representative magnitude characterizing this profile, and 
obtaining the static pressure using the prescribed total pressure. The 
guessed representative velocity magnitude, U a , is then updated as the 

overall solution evolves. For channel flows, McDonald and Briley [27] 
have suggested updating the total pressure distribution within the 
boundary layer, in order to maintain the required velocity- and 
temperature-profile shapes prescribed at the inflow boundary. This 
implies that for fully developed flow conditions at the inlet, as used 
in the present channel-flow studies, the total pressure distribution has 
to be updated over the entire channel width. 

The procedure for updating the velocity profile for cascade flows 
will be described in the next chapter. 

2.4.2 Wall-Wall Boundary Conditions 

The governing equations given by Eq. (2.18) consist of one first- 
order equation, namely, the continuity equation, and three second-order 
equations, namely, the x- and y-momentum and the energy equations. 
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Here, the order refers to the highest-order derivative in the n 
direction. . These comprise a system of seventh order with respect to n, 
so that a total of seven boundary conditions need to be specified along 
the two n ■ constant boundaries. At the wall, the no-slip boundary 
condition is imposed; also, the walls are assumed to be impermeable and, 
hence, there is no injection or suction at the surface. In addition, 
the temperature, T w , at the surface is specified. These constitute a 

total of six boundary conditions at the two surfaces. Therefore, one 
additional condition has yet to be specified. A valid flow approxima- 
tion such as Op/3n) » 0 can be imposed as an additional boundary 
condition. The resulting solution will reflect the approximations 
inherent in the boundary condition. Another method to obtain the 
additional condition is to write the governing equations in one-sided 
difference form at the wall, as has been done by Rubin and Lin [29] and 
Briley and McDonald [30]. 

In the present study, an approximate form of the normal momentum 
equation, obtained by dropping the viscous terms in that equation and 
written at the first cell center near the wall surface, is used as the 
additional condition. The viscous terms in the normal momentum equation 
can be shown to be negligible near the wall surface for most of the 
flows considered in this research. 

In order to ascertain that enough independent equations are 
available at a particular streamwise location, a typical grid line along 
the n direction, consisting of five computational points as shown in 
Fig. 2a, is considered. Counting four unknowns, namely, 
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T 

Q - (p, pu, pv, pe ) , per point, this line involves a total of 20 

v 

unknowns. The independent equations that can be written are the 
continuity equation written at the 4 points denoted by c in Fig. 2a, and 
the momentum equations and the energy equation written at the 3 points 
denoted by x, y and e, respectively. This constitutes a total of 13 
independent equations. When combined with the six specified wall- 
boundary conditions, these equations leave one additional condition to 
be specified. This additional condition is taken to be the reduced y- 
momentum equation as discussed in this section, and may be imposed near 
either wall surface. 

These boundary conditions can be stated mathematically as follows. 

At the walls, 

u = u =0 , 

w 

v = v =0 
w 

and T = T . (2 - 21 5 

w 

The reduced normal momentum equation near one of the walls can be 
written as 

j (pvU + 5 y P)} + { j (P vV + n y P)1 = 0 . (2.22) 

2.4.3 Wake-wake Boundary Conditions 

For cascade flows, specification of boundary conditions along the 
wake boundary needs to be considered. For symmetric configurations, the 
wake boundary, is in addition to being the wake centerline, also a line 
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of symmetry. For unstaggered cascades, a 'line-periodic' grid is 
employed. In this type of grid, the same n - coordinate line connects the 
corresponding periodic points, such as points 1 and 5 or 0 and 4 in 
Fig. 2b. In this case, the periodic-boundary conditions can be enforced 
implicitly. The periodicity condition is that the corresponding values 

T 

of all the flow variables, Q - (p, pu, pv, pe fc ) > and the normal 

derivatives, u , v and T , of the velocities and temperature which are 
n n n 

governed by second-order equations, namely, the momentum and energy 
equations, must be the same at corresponding periodic points along the 
wake boundaries. It should be mentioned that, in terms of the conserved 

T 

variables, Q * (p, pu, pv, pe t ) , the repeating condition on the n- 

derivatives must be satisfied for all four elements of Q . This 

n 

condition can be written as 

Qi - Q s or Q 0 = Q„ (2.23) 

for a typical computational line consisting of points 1 through 5 as 
shown in Fig. 2b. Imposing the periodicity boundary condition described 
above between points 1 and 5, leaves a total of 16 unknowns counting 
four unknown variables per computational point. The system of equations 
that can be written along this computational line consists of the 
continuity equation at 5 points, and the momentum and energy equations 
at 4 points, denoted by c, x, y and e, respectively, in Fig. 2b. As 
periodicity has already been imposed, the continuity equation, which is 
of the first order, written employing points 4 and 5, becomes identical 
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to that employing points 0 and 1 . Hence, there are only 16 independent 
equations and the system is closed. Also, in the case where symmetry 
exists along the wake boundary, the viscous terms in the y-momentum 
equation can be dropped, and this reduced first-order equation can be 
written employing points 1 and 2. 

2.4.4 Wall-Wake Boundary Conditions 

This type of boundary condition is needed for cascades with stagger 
when a 'region-periodic' grid is employed. A 'region-periodic' grid is 
one in which the corresponding periodic points in the flow are not 
connected by the same n-coordinate line. This type of grid has to be 
employed for cascades with large stagger in order to avoid excessive 
skewness of the coordinates. The use of this type of coordinates, in 
conjunction with a marching procedure, forces the periodicity condition 
to be imposed in an explicit manner. 

Figure 2c shows a typical grid consisting of six points along an n - 
coordinate line in the wall-wake region. The 18 independent equations 
along this line consist of the continuity equation written at 5 points 
and the momentum and energy equations written at 4 points, in addition 
to the reduced momentum equation written at the wall surface. The six 
boundary conditions consist of the zero slip, zero inject ion/ suet ion and 
the specified temperature at the wall surface and the specified velocity 
and temperature conditions at the wake boundary. The conditions at the 
wake boundary are obtained from the corresponding periodic point along 
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the upper wake boundary in the flow. These conditions can be stated 
mathematically as follows. 

At the wall, the conditions are 
u - u w - 0 , 

v - v w - 0 

and T - T (2.24) 

w 

Along the wake, for example, at point 0 in Fig. 2c, the conditions are 

u ° “ u a ' 

Vo - v a 

and T 0 - T (2.25) 

u a 

Here, 'O' and 'a' are the corresponding periodic points. 

2.4.5 Wake-wake Boundary Conditions (Region-Periodic Grid) 

Along this type of boundary for staggered cascades employing a 
region-periodic grid, the periodicity condition is imposed in an 
explicit manner. The independent equations to be considered along a 
typical computational line consisting of six grid points 0 through 5 are 
the continuity equation written at 5 points and the (x,y) momentum and 
energy equations written at 4 points. The boundary conditions consist 

T 

of the specified values of Q - (p, pu, pv, pe fc ) at point 5, obtained 

from the values at the corresponding periodic point in the flow, such as 

point b in Fig. 2d. For this purpose, the most recent values of Q are 
used. This implies that only three more boundary conditions can be 
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supplied for the seventh-order system. Hence, the velocities u and v 
and the temperature, which are lagged in time, are specified at point 0. 
These conditions can be stated as follows. 


Along the wake, at point 5, 


P, - P b • 


u. - u v 


v g - v b and T. 


and at point 0, 


u o - U a 


and T 0 ■ 


(2.26) 


(2.27) 
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CHAPTER 3 


NUMERICAL PROCEDURE 

The Interacting parabolized Navier-Stokes equations (2.18), 
described in the preceding chapter, are a set of nonlinear coupled 
partial differential equations. Analytical solution of this system of 
equations exists only for a small, special class of problems. Hence, 
for the general problems of present interest, a numerical solution of 
these equations is sought. 

The linearization and discretization of the governing equations, 
the solution procedure for the discretized set of equations, the 
implementation of the boundary conditions, including the periodicity 
boundary condition f.r cascade flows, as well as the treatment of 
problems with flow separation are detailed in the following sections of 
this chapter. 


3-1 Quasi-linearization and Discretization 

The system of governing equations (2.18) is first re-written here 
for easy reference. 


3E 3F 



(3.1 ) 


where 



and F 



(3.2) 
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As mentioned above, this system of equations is nonlinear. 

Therefore, these equations must be linearized or quasi-linearized in 
order to obtain a system of algebraic equations amenable to numerical 
solution. In the present study, employing forward marching for the 

solution vector Q, a quasi-linearization of the nonlinear terms at a 
given streamwise station is carried out about the solution at the 
preceding streamwise location. From Eq. (3.2), it is seen that the 
quasi-linearization of the inviscid flux vectors E and F requires that 

the flux vectors E and F be quasi-linearized. This is achieved by using 


Taylor’s series expansions. The results can be expressed as 

follows. 

E i+1 - E 1 + 

(3.3) 

F i+1 - F 1 + B 1 A*Q 

(3 • *0 

where 


4*0 - Q 1 * 1 - Q 1 . 


J 1 . il‘ and 5 1 - 4 1 • 

(3-5 a-c) 

3Q 30 



Here, the superscripts i and (i+1) denote two successive streamwise 
locations as shown in Fig. 3- The quantities with superscript (i*1) are 
the unknown terms, which contribute to the nonlinearities in the 


equations. 

Recognizing the fact that the inviscid flux vectors, E (Q) and 

F (Q) are homogeneous functions of degree 1 In Q, and using the property 
of homogeneous functions (see Ref. [31]), one can write 
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“ Q ■ E and Q - F 

3Q 9q (3.6) 

This property allows some conceptual simplifications of the equations 

and leads to computational efficiency. Using this property of 

homogeneous functions, Eqs. (3.3) and (3.H) can be written as follows. 

=i + 1 ri -i+1 
E - A Q , 

and F 1 * 1 - B 1 Q i+1 

* (3.7 a,b) 

The Jacobian matrices A and B, defined in Eqs. ( 3 .5 b.c). are given 
as follows. 


(3E/3Q) = 


4> 2 - u 2 


-uv 


u ( 2<f> 2 -Ye t ) 


and 


1 ' 0 1 0 
(3-Y)u ; “(Y-1 ) v ' (Y-1 ) 

I I 

V u 0 

i i 

(7et _ 4> 2 j -(Y-1 )uv 1 Yu 
-(Y-1)u 2 }| 


B - (3F/3Q) 


0 

0 

' 1 

| 

1 0 

-uv 

V 

~l ~ u ~ ’ 
1 

i 0 

l 

<t > 2 - u 2 

-(Y-1 )u 

1 (3~y)v 

i _ 

1 (Y-1 ) 

_ 

(2<f> 2 -ye t ) 

-(Y-1 )uv 

r (T v»‘ 

i 

1 Yu 
1 

t 


i 

1 -(Y-1 )u 2 } 

1 

1 



j 

— 
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The Jacobian matrix A 1 contains metric quantities evaluated at station 
i. The - on the Jacobian matrices A and B in Eq. (3.9) is used to 
denote that the metrics in these matrices are evaluated at station 
(i+1), and distinguishes these matrices from the Jacobian matrix in 
Eq. (3.10). 
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The viscous term F y is quasi-linearized as follows. A typical term 

38 

of F is of the form (a ■r 2 ’). Starting with the Taylor's series 

V oH 

38 i + 1 

expansion for (a — ) and re-arranging leads to the expression 

(a|^) i+1 - a 1+1 (l^) 1 + a i+1 (96/3Q A*Q) (3.11 a) 

3n Bn 

Further, it is observed that the term 8 is homogeneous of degree zero, 
in Q. Hence, 


(38/3Q) Q - 0. 

Therefore, Eq. (3.11 a) becomes 

, 36 v i + 1 i + 1 / 36 s i ^ i + 1 n i + 1 

(a = a (t^) + a (36/3Q) Q 


(3.11 b) 


(3.11 c) 


Using Eq. (3.11 c), the viscous term F^ +1 can be written as 
F 1 *' - . i 1 Q 1 *' 

V V V 

where a typical term of F* is of the form a* ^ (-Ir^)^ , and of ^ is 

1V v 3n v 

a i+1 (36/3Q) i Q i+1 • Again, the - denotes that the metrics involved 
are evaluated at station (i+1). 

The discretization of the derivative terms is discussed next. The 

3 E i + 1 

^-derivative term (-r— ) is considered first. This term can be 


written as 
i+1 


f 3E^ ' a r3 £ i i+1 x ^ f 9E i i 

^3f/ 01 ^35^ (1 01 ) [g^J 


( 3 . 12 ) 
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The values of 9i are chosen depending on whether an implicit, explicit 
or the Crank-Nicolson scheme is to be employed. The various values of 
commonly employed, and the corresponding schemes, are 

( 1 for an implicit scheme , 

0 for an explicit scheme , 

1/2 for Crank-Nicolson scheme . (3-13) 

Also, depending on the type of differencing and the desired 
accuracy, one can write 


ap 1+1 A 1 E A 1_1 E 

(||) - o*e 2 ) V - «. -sr- 


(3.1*0 


where A 1 is the forward difference operator defined by Eq. (3-5 a). 

The type of differencing, and the corresponding order of accuracy, for 

various values of 0 2 are given as 

{ -1/2, central, 0(AE, 2 ) 

0 , two-point backward, 0(A£) 

1/2 , three-point backward, 0(AE 2 ) (3-15 

Combining Eqs. (3.12) and (3.1*1) leads to, 


0i (|f) i+1 * O-0i ) (ff) 1 - 0*e.> 4r * 02 TT 


(3-16) 


An explicit scheme usually has some stability condition such as the 
Courant-Friedrichs-Lewy (CFL) condition associated with it. Hence, an 
implicit scheme, corresponding to 0j - 1, is employed in the present 
work, to enhance the numerical stability of the resulting method. To 
facilitate a marching type of procedure, either two-point backward 
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differencing or three-point backward differencing has to be employed. 

In the present work, the two-point backward differencing, with e 2 - 0, 
has been employed. 

3E 

Substituting for from the governing differential equation 

(3.1) into Eq. (3.16) gives 


" 6l h (F “ F v )i + 1 + < 1 - 8 |> (|f) - ( 1 + 0a ) - 0, A 1 ’ 1 E 


AC 


AC 


.(3.17 a) 


Equation (3.17 a) can be rearranged as follows: 


A 1 E 


TTt 4 1 ' 1 E - rrt; « Is < F - F /*’ 


- nfi 45 h <F - F v >1 


(3.17 b) 


Substituting for the quantities at station (i+1) using the quasi- 
linearization given by Eqs. (3.9) and (3.11), Eq. ( 3.17 b) can be 
rewritten as 


A 1 Q 1+1 - E 1 


E -T?t {(§i ‘ § v 1} " F v i} 


- 45 rift h <F - F v )1 


(3.18) 


Equation (3.18) can be re-arra.nged to contain the unknowns at station 
(i+1) in the incremental or delta form, by subtracting 

li ~ q1 + fe {(§i - 5 1 ) 

from both sides of Eq. (3.18) and is given as 
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A 1 * « TTaT In l(§1 ' 4151 


- E 1 - A 1 Q 1 + -rri 


i-1 


1-0, 


A 1 Q 1 * tHt a B - A5 TTft fe (F ‘V 


9 1 3 \ 

T^eT ^ {( 


B * - Bj) Q‘ - F*} 


Here, 


(3.19 a) 


A 1 Q 1 


B 1 Q 1 


(/) 


x^i+l si 


*- ♦ (/) 


1+1 F 1 


(^) i + 1 I 1 - (^) l + 1 F 1 


and B* Q L 


0. (3.19 t-d) 

It may be important to recall that Eqs. (3*19 b,c) are obtained 


from the definitions of the Jacobian matrices A and B, given by Eqs, 
(3.9 b,d), and by using the homogeneous property [Eqs. (3-6)] of the 


flux vectors E and F. 

For an implicit scheme employing a two-point backward difference, 

-1 and e 2 - 0. Hence, Eq. (3-19 a) reduces to the following. 

A 1 A*Q ♦ AS CCB 1 “B*) A X Q} 

- E 1 - A 1 Q 1 - AC (B 1 Q 1 ~ fJ) (3-20) 

The next step in the numerical procedure is the introduction of 
upstream influence. For flows with strong viscous-inviscid interaction, 
as mentioned in Section 2.3, the streamwise pressure gradient term is 
split according to Eq. (2.20) and upstream influence is introduced by 
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forward differencing the ’(1-w) 3p/3E * term. With this *« - split’, 
the governing equation (3.1) can be written as 


3E 


3F 

3n 




(3.21 a) 


where 


a 



0 

? x,f 
(1-») C y>f 


L 


0 


(3.21 b) 


The subscript f on the pressure gradient term denotes that this 

term is forward differenced and those on the metrics E and E denote 

x y 

that these are obtained using coordinate values at station (i+2). This 
representation of the metric terms was arrived at by applying the 
procedure for a test case in which a fully developed flow was reproduced 
in a straight channel using a coordinate system with metrics varying 


along the streamwise direction. 

* 

The flux vector E is different from E given in Eq. (2.16), in that 
the pressure terms in it are multiplied by the factor u>. It is given as 

pU *] 

* 1 

E - — puU + E x up j 


pvU + Ey up 
(pe fc +p) U 


( 3 . 22 ) 


where U is the contravariant velocity given by Eq. (2.17 a). 
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For low subsonic flows considered in the present study, according 
to Eq. (2.19), u - 0. This implies that the total streamwise pressure- 
gradient term is forward differenced. Also, it is clear from 
Eq. (3.21 a) that the streamwise pressure gradient term is no longer in 
the SCL form. This necessitates that the pressure gradient term p^ 


should no longer be in the SCL form, as the starting governing equations 
were written in the conservation form in Cartesian coordinates. 
Accordingly, Eq. (3.21 a) must be written in the following form. 


3E* , 3F * h 3 P - i (f ) 

H aS aS an » 


- a 


12 


(3-23 a) 


where 


b 



0 


( 1 -<d) n x 
(i-id) n 

i 

0 I 

-i 


(3-23 b) 


and F* - F - b p . O 

The procedure for quasi-linearization of Eq. (3*23 a) is similar to 
that used to obtain Eq. (3-20). The resulting equation has an extra 
source term on the right hand side, arising out of the pressure gradient 
(3p/3£)f term. The quasi-linearized form of Eq. (3.23 a) becomes 

A * 1 * 4 ‘Q . M [ & US' 1 -Sj> A l Q ) * b 1 *’ f; «%* A*Q] 


E* 1 - A* 1 Q 1 - A5 (f^ (B* 1 0 1 ) * b 1 * 1 (I?) 1 - I: (fi)} 


3n 


3n v' 


i + 1 A i+1 n 

- a A P 


(3.24) 
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The asterisk on the Jacobian matrices A and B denotes that these are 

ft ft 

obtained from the corresponding flux vectors E and F , respectively. 

A possible method of solution of Eq. (3.24) is to employ an 
alternating-direction explicit (ADE) method. In the first step of this 

method, A Q is computed starting at the inflow boundary and proceeding 

towards the outflow boundary and in the second step, A i+1 p is computed 
starting from the outflow boundary and proceeding towards the inflow 
boundary. Such a method requires the inclusion of a time— derivative 
term, p fc , in the momentum equations containing 0p/30 f term, in order 


to unlock the solution from its initial conditions. This method, like 
the pressure updating procedure used previously by K. Ghia and U. Ghia 
[33], is capable of transmitting the downstream disturbances upstream 
efficiently. It is, however, algebraically much simpler and has been 
successfully used by Barnett and Davis [34], for solving supersonic 
external flow problems. 

With the inclusion of the time-derivative term, the governing 
equation (3.23 a) assumes the following form. 

* » 


3E 3F 

35 3n 


♦ b 


l£ _ 3_ 

3n 3n 


<v 


“ a 


IE 

35, 


+ a 


3£ 

3t 


(3-25) 


The quasilinearized form of Eq. (3.25) is obtained by utilizing Eq 
(3.24), in addition to discretizing the (3p/3t) term as 

an 1 + 1 

(|f) 


1 


At 


(p i+1 ’ n+1/2 - p i+1 ’ n ) 


(3.26) 
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, , i+l.n+1/2 . . i.n+1/2 .... 

and quasi linearizing p about p . The resulting 

equation is 

A* 1 A*Q ♦ A*- [ ~ {(B* 1 - B^) A*Q} + b 1 + 1 t(^) 1 A^} 


an 


3Q 


ai + 1 {(&) 1 A*Q }] 

aQ 


At 


*i -*i -i .3 .~*i -i . . i + 1 3p 3 ,=ix, 

- E 1 - A 1 « - AC (B Q ) * b ^ - 5S <F y » 


. a i + 1 ,M (p i."* 1/2 - p 1 *'-") - 4 ‘*' p") 
At 


(3-27) 


The discretization of the above equation in the n direction is 
considered next. The discretization is performed such that the 
discretized form of the inviscid portion of the equations constitute a 
consistent set of equations by themselves; the viscous portion of the 
equation is also discretized in a self-consistent manner. The inviscid 
part of Eq. (3-27) is given by 

A* 1 A X Q + AC [ UB* 1 A 1 ?} + b 1 + 1 A X Q 


3n 


3Q 


- a1q} ] 
At aq 


*i ~*i -i ,3 .-»i riv .i + 1 3p , 

■ E 1 - A Q - A5 <B q ) * b ^ ) 

. a 1 ” {if (p 1 ’"* 172 - p 1 * 1 '") - A U ' P) • 


( 3 . 28 ) 


All normal derivatives are represented by central differences, with 
the first - order derivative representation involving points across one 
normal mesh interval only. The discretized equations are second-order 


35 



accurate in the n direction. Equation (3.28) is written at a normal 
mesh midpoint, j + 1/2 , as 


\ A*Q )j + (A * 1 A*Q ) J+1 J ♦ [( B * 1 A*Q ) j+1 - (B * A a Q). 


i I- 


* -S'* "I - <1 


' rh vU V, * ‘I 1 * 15 V 

i I Bj *‘ * E*;, - ( A * 1 Q 1 ^ - ( A * 1 Q 1 ) } 


- { ( B * 1 Q 1 )... " (B* 1 Q 1 ), + b 1 + 1 fn 1 -" - n 1 ’ 11 

An j + i 


J -J../2 - »j’“» 


1 o i+1 r r / i,n+1/2 i + 1 ,n, 

2 a j+1 /2 ^ IT {( P - P > 


j + T 


+ (P 


i,n+1/2 i + 1 ,n,. , .i + 1 n ,i+1 ni 
- P )j) - A Pj„ - A Pj] 


(3.29 a) 


where 


i+1 1 , i+1 i+1 . 

a j.,/2-2 U J * a j.1> 


and 


b 


i + 1 

j + 1/2 


2 



b i+ !) 

j + i 


(3.29 b.c) 


In order to include the viscous terms at mesh point j, the inviscid 
part of the equation should also be written at j. This is obtained by 
writing a discretized equation similar to Eq. (3*29) at location (j-1/2) 
and forming the arithmetic mean of this equation and Eq. (3.29). The 
resulting equation is 


36 


(A* 1 A^Q ) j+1 + l (A* 1 A*Q )j + { ( A* 1 A l Q) ^ 


AC 

2 An 


[(B* 1 A*Q ) j+1 - (B A A x Q) j _ 1 


**i <- 


‘rft t a ]*i/2 1( ^ 415 >j.i * ‘‘o V 


+ b j + 1/2 {( ^ A Q } j+1 


- (3H a a q ).} 

3Q J 

i 


*>}I] /2 {(% A 1 Q) - (^£ A 1 Q) }] 

J 3Q J 9Q J 1 


9Q 

i 


+ a i + 1/2 l( ^ a1q) 1 + ( ^ i A 1 Q) 1 .}] 

J 3Q J 3Q J 

I 1 E*i, '(A * 1 Q 1 )j +1 * E*i, -(A* 1 Q^j.,} * | (Ej* 1 -(A* 1 Q 1 )j > 

-rh f< 6,1 5i) j*i - (5n 


i+1 , i,n i.n, 

Vl/2 ( Pj .1 “ »J 


i + 1 


i.n i.n, 


. 4. / Jl 9 i i i • 1 1 \ . 1 ' 1 / l • ii 1 • ii v > 

* b - (P <- - P. > * Pj-,/2 (P j - Pj-, )J 


. 1 i+1 AC ,, i,n+1/2 i + 1,ru , i t n+1/2 i+1 ,n, , 

* « a j+i /z it l(p - p >j.i * <p - p y 

_ 1 i+1 AC ,, i,n+1/2 i+1 ,n, , i,n+1/2 i+1,n, , 

* 1 a j-i/2 it {( p - p >j * ( p - p y,> 


1 r i+1 ,, 1+1 n . ,i+i n. 1+1 , A i+i n . 1+1 n . , 

*t 1 a j + 1/2 4 PjH-1 * 4 Pj> * a J-1/2 (4 Pj* 4 Pj-y 


i + 1 n 


i+1 n 


i + 1 


i+1 n .i+1 n 


(3.30) 
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The discretization of the viscous terms in Eq. (3.27), for example, 
3(B* A*Q)/3n , is considered next. The viscous flux vector F y and the 


Jacobian matrix B y can be given as follows: 


a, u + a, v 
1 n 3 n 


J Re 


a, u + a, v 


a, uu + a, w + a, (uv) 

1 n n n 


♦ l., (e. - 


u 2 + v 2 


(3.31 a 




- a 3 (-) 

p n 

- a 3 (-) 

3 p n 


a» (-) 

1 p n 


1 a 3 (1) 

3 p n 


! o 
'•o 


where 


I 


J Re ! 


- a 3 (J) 

p n 

- a 2 (-) 

2 p n 


a 3 (~) 

p n 


a 2 (t) 


p n 




- a 3 (— ) 

1 p n 

V 2 

" <7-> n 

p n 

- 2a 3 (^) 

p n 


‘ ^ ( ?>n 


(a,-aj (J) n 
+a, (J) 

p n 


a 3 (^) ! 

p n 

*<».-*.) <?>„ <?> 


e - e t - (u 2 +v 2 ) , 


*• ■ | V * V 


38 



2 


2 


** ' * 5 \ 


l« ■ <% n y )/3 

and 




Y_ 

Pr 



(3-31 b-g) 


It should be mentioned that when the term (B y AQ) is formed, the AQ terra 

should be contained in the n~deri vatives appearing in Eq. (3*31 b). A 

3 ( 86 ) 

typical term of this product can be expressed as [a . For 

example, for the element corresponding to the (2,2) location of 
( B v AQ), 0=1, , B = and <5 - A(pv). The discretization of a 

typical term — {a -r- (B 6)} is performed by evaluating the quantity 
an on 

in brackets at two successive mesh mid-points, such as (j+1/2) and 
(j-1/2), and forming a difference expression at the mesh point j to 
obtain a second-order accurate representation for this term. This is 
outlined below. 


3 _ 

3n 


(o (B 6)} . 

3n j 


An 




(B 6)) 


j + 1/2 


(a 9n (B $))j_i /2 } 


where 


and 


(a 3fi (6 6 ) ) j + i /2 " l a j + 1/2 (B j + i 6 j + 1 " B j 6 j^ 


(a 3n (6 «))j_ 1/2 “ An { a j--i/ 2 (6 j 6 j " B j-1 6 j-1^ * 


(3.32 a-c) 
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Using the discretization described by Eqs. (3-30) and (3.32), Eq. 
(3.27) can be written in a compact form as 

A J 4 ‘V, * B J “S * C J 4 ‘V, ’ D j • (3-33) 

The discretization described above applies at a general interior 
point; the treatment of boundary points will be discussed in the next 
section. 

Equation (3.33), written at all the mesh locations j along a line 

(i+1), results in a block-tridiagonal system of equations, with A., B 

<3 J 

and Cj being (4x4) matrices. This system of equations, can then be 

solved using L-U decomposition of the coefficient matrix of the system. 
The implicit solution procedure for such a 3y3tem can be found in 
Ref. [35], by Anderson, Tannehill and Pletcher. 

3.2 Details of the Solution Procedure 

The solution procedure consists of two time steps. In the first 
step, the solution proceeds from the inflow boundary towards the outflow 

boundary, employing Eq. (3.25) to update Q . The pressure field is 
updated in the second step of the procedure. This step proceeds from 
the outflow boundary towards the inflow boundary. The outflow boundary 
condition on pressure is directly employed during this step. These two 
steps can be expressed as follows. 

Step,: R n *> /2 . - a (f|)J ♦ a (|f )"* ,/2 (3.3- a, 
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o „n+ 1/2 cSp^n+l r3p-.n+1 

Step 2: R - - a [^|J f + a [-^J ( 3.34 b ) 

where R contains all terms in Eq. (3.25) except the pressure terms which 
appear explicitly in the above equation. In Eq. (3-34), superscripts 
denote time levels. Thus, for example, the superscript (n+1/2) denotes 
that the time derivative is evaluated at time level (n+1/2). Backward- 
difference approximations are used for the time derivatives, so that the 
time derivative at (n+ 1 / 2 ) employs the pressure at time levels n and 
(n+ 1 / 2 ). 


A simpler equation for the second step can be obtained by 

eliminating R 0 * 1 ^ 2 between equations (3-3*1 a) and (3-3*1 b). The 

resulting equation is 

flE) n + |'9p-in + 1/2 _ _ j-a^n +1 , 3p^n+1 

l 3C J f 3t l 3C J f 3t 

The discretized form of Eq. (3-35) at station (i+1) is 

, i+ 2 ,n i+ 1 ,n, AC . i+ 1 , n+ 1/2 i+ 1 ,n. 

- (P " P ) + tt(p - P ) 


(3.35) 


, i+2,n+1 i + 1,n+1. AC , i + 1 ,n+1 
(p ’ - P ) + S <P 


i+ 1 , n+ 1/2 
P 


(3-36) 


From this equation, 


i+ 1 ,n +1 
P 


can be solved for in an explicit manner as 


pi+ 1 ,n +1 _ [ p i+ 1 ,n _ p i + 2 ,n + p i+ 2 ,n +1 

+ || ( 2 p i+1,n+1/2 - p i+1,n )] / (l+ ||) (3-37) 

Equation (3-37) is applied along lines of constant n and the 
pressure field is updated by marching upstream. The prescribed 
condition on the pressure at the outflow boundary is imposed via the 
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pi+2,n+1 tepm (3.37). The two steps described above constitute 

one global iteration, during which the flow solution is advanced from 
the time level n to time level (n+1). 

If the pressure terms are retained in the SCL form, then the 
equation for the upstream marching step, corresponding to Eq. (3.35), 
will be a set of two equations in the single variable p. This system 
can be reduced to an equation similar to Eq. (3.35) by combining the two 
equations after multiplying each of them by the appropriate metrics. 

This is equivalent to taking a projection of the two equations along a 
line of constant n. 


3 . 3 Updating of Velocity Profile at the Inflow Boundary 

The reference velocity, U ay g, characterizing the velocity profile 

at the inlet, needs to be updated before the next global iteratirn is 
performed. This is necessary because at the inflow boundary the total 
pressure and temperature are prescribed as boundary conditions, so that 
U av g has to be guessed to initiate the solution procedure. The updating 

of U av g is as follows. At the inflow boundary, knowing the dimensional 

* * 

total pressure p^ and the dimensional static pressure p (which is 

evaluated through the upstream marching step), the local Mach number M 
can be obtained using the isentropic relation 


p 


(3.38 a) 
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Then, using the dimensional total temperature T Q prescribed at the 

inflow boundary, the static temperature T can be obtained from the 
following relation 

i - ■ (3 ' 38M 

T 

The dimensional density p* can then be found, using the equation of 
state, as 

# * * (3-38 c) 

p - p /R T 

Knowing M* and T*. together with the definition of the local sonio 

# *2 *2 1 / 2 

velocity, the local dimensional speed V (■ tu + v } ) , can be 


determined from the relation 


* , # 
/ 2 = M 2 YRT 


(3.38 d) 


* * 

The Cartesian components of the velocity u and v can e 
determined from the given flow direction, 0, at the inlet so that 

u* . iv* /(l.tan*0»' /2 < 3 ' 38 e> 

and 


* * (3.38 f) 

v - u tan0 

The reference velocity U ayg is then obtained using Eq. (2.6) and 


all the variables are then re-nondimens ionalized. 
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3.4 Convergence Criteria 


The two steps of the solution procedure described in Section 3 2 
are repeated until convergence is achieved. To test for convergence, 
the maximum absolute value of the error (L^-norm) and the root-mean 


square of the relative error (L 2 -norm) in the pressure field are 
monitored. These are defined, respectively, as 


'abs 


and 


max 

- i-I.IMAX 
j-1 , JMAX 

JMAX IMAX 


p j."- . Pj ,n 


JMAX IMAX p* ,n 

rel " ( ^ l ( 1_ i n+1 ) 2 } 1/2 /(At IMAX JMAX) (3.39 a,b! 

j-1 i-1 P ’ 


Convergence is said to have been achieved when e < io ^ and 

abs 


e rel S 10 


-6 


3 * 5 Imple mentation of the Boundary Conditions along n J Constant 
Boundaries 

In the present study, the boundary conditions at the walls and, in 
the case of cascade flows, along the wake boundaries are implemented in 
an implicit manner, consistent with the numerical procedure employed in 
the interior of the computational domain. The implicit treatment of 
these boundary conditions in an otherwise already implicit solution 
procedure removes the mesh spacing constraints encountered in an 
explicit scheme and also aids in enhancing the convergence process. 
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3.5.1 Wall Boundary Condition 


As mentioned in Section 2.H.2, the zero-slip and zero suction/ 
injection conditions at the wall, together with the wall temperature, 
are specified. These are written as 
u - u w - 0, 

v » v - 0 
w 

and T - T <3- HO a-c) 

w 


Expressed in terms of the increments in the variable 


Q - 

(p. pu, 

pv» 

T 

pe fc ) , Eq. 

(3. HO) yields 



i + 1 

- u 

A*P 

+ A 1 ( pu) = 

i . i 

p A u , 



i+1 

- V 

A X P 

+ A 1 ( p v) = 

i . i 
p A v 


and 

i+1 
' e t 


+ A 1 (pe t ) 

i . i 

“ P A e t * 

(3. HI a-c) 


In the above equation, all quantities with superscript (i+D are 
known from the conditions given by Eq. (3. HO) and the right-hand side of 


Eq. (3. HI) can be evaluated using the known solution vector Q 1 at 
station i. The three equations given by Eq. (3. Hi), together with 
either the continuity equation or the reduced y-momentum equation, 
constitute the four equations at the wall boundaries. 

3.5.2 Periodicity Boundary Condition for Cascades 

For flows through cascades employing a 'line-periodic* grid, the 
periodicity boundary condition can be imposed in an implicit manner. 

The periodicity condition, as described in Section 2.H.3, requires that 
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the flow variables have the same values at corresponding periodic points 
along the wake boundaries. When this condition is imposed, the 

discretized equation (3-33) assumes the following form. 

At j»1, 

A 1 * B 1 ' C , ‘ D , «•«) 
and, at j-JMAX-1 , 


JMAX-1 




JMAX-2 


+ B 


JMAX-1 


A 1 Q 


JMAX-1 


+ C JMAX-1 A Q 1 


JMAX-1 


(3.^3) 


Hore, JMAX is the index corresponding to the maximum value of n. 
Equations (3.42) and (3.43), together with Eq. (3-33) written at each 
interior normal mesh point j=2 through JMAX-2, form a system of 
equations which is basically a tridiagonal system, except for non-zero 
corner elements. The corresponding coefficient matrix is shown below. 








r 



B 1 

C 1 

0 

0 

A 1 


AO, 


D 1 

A 2 

B 2 

C 2 

0 

0 


aq 2 


°2 

• 

• 

• 

• 

• 


• 


• 

* 

• 

• 

• 

• 


1 

- 

» 

0 

0 

0 

A JMAX-2 B JMAX-2 

C JMAX-2 


A ® JMAX-2 


°JMAX-2 

C JMAX-1 

0 


A JMAX-1 

B JMAX-1 


A ^ JMAX-1 ■ 


° JMAX-1 







* — 


... 


(3.44) 
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Here, A, B and C are (4x4) blocks, while AQ annd D are (4x1) blocks. 

The procedure for solving the above periodic block tridiagonal system in 
an implicit (non-iterative) manner is detailed in Appendix A. 

3.6 Separated Flow Modeling 

In this section, the approximations involved in obtaining the 
governing equations in separated flow regions are described. It is 
known that forward marching in space with the parabolized Navier-Stokes 
equations in regions of reversed flow, that is, where the tangential 
contravariant velocity component U is negative, is unstable. This 
instability can be overcome if all the equations are forward 
differenced in the regions of reverse flow. This requires that, in 
addition to the pressure, the preceding iterate of the solution vector 

Q be stored in these regions . 

Reyhner and Flugge-Lotz [ 36 ] have suggested a simple alternative to 
this situation. They suggested that, in the reverse-flow region, the 
convective term u 3u/3x in the momentum equation be represented by C |u| 
3u/3x, where C is zero or a small positive constant. This 
representation, known as the FLARE approximation, assumes that the 
convective terms are small in regions of reverse flow and is valid when 
the reverse flow velocities are small. This approximation is employed 
in the present study also, by neglecting all the convective terms in the 
momentum and energy equations in the reverse-flow region. Hence, in 
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regions of reverse flow, the governing equations take the following 


form. 

When U - (£ u + 5 v)<0 , 

x y 

35 3n 3ri ' »' 

where 


pu 


pV 

5 x P 

' F ~1 

n p 
x K 

5 V P 


n„ p 

y 


y 

0 


0 


and F v is as described earlier in Eq. (2.15). 


(3.^5) 


3-7 Discretization of Metric Terms 

The numerical representation of the metric coefficients arising 

due to a general coordinate transformation from the physical domain to a 

computational domain, given by Eq. (2.10), is described in this section. 

The metric coefficients, such as E , E etc., are obtained from the 

*x *y 

derivatives x , y , etc., using the relation given by Eq. (2.13). 

n n 

The discretization of the metric derivatives should be done in a 
manner consistent with the discretization of the governing equations. 
Hindman [37] has shown the appropriate representation of the metrics for 
solving a 1-dimensional wave equation, using MacCormack's scheme and 
various forms of the governing equation, such as the strong- 
conservation-law (SCL) form, the weak-conservat ion-law (WCL) form, etc. 
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In his work, he employed a simple test of reproducing uniform flow, 
starting with the entire computational mesh Initialized with a uniform 
flou and advancing In time, employing a selected numerical Integration 
algorithm. In the present study a similar test Is performed to 
determine the appropriate representation of the metric derivatives and 
ls detailed In Appendix B. As shown In the Appendix B. representation 
Of the transformed eduations in the SCL form retires the following 

relation 

E * n v 

!_(_£)♦ (-p) - o 

35 '■J ' 3n J 


and 


8_ (!i) . f & - 0 . 

35 ' J 3n J 


he satisfied In the discretized form. This Implies that the discretized 

representation of the n derivative In ““ “ 

3 \. also, the ^-derivatives 

that employed for the n-derivative in 3n (y^). 

o s are summarized 

are to he discretized In a similar manner. The results 
here . 

. . (i ax the coordinate derivatives take the 
At an interior point (i*j)> 

following form. 


x n|i,J ' <X i.j*i ' 

y n|l.J ‘ 


)/2An 


)/2An , 
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X 


i.J 


and 


- V,J> 


X l-T , j-1 


)] 


■ *ac * ( J r t. j-i - 

, 1 ( 

2A5 y i,j " y i-i • (3.46 a _ d) 

From Eq. (3.H6), It is clear that the n-derivatives are repreaented 
by second-order aocurate central differences and the {-derivatives are 
represented as averages of first-order accurate backward differences. 

At a boundary point J-JMIN or JMAX, the {-derivatives are 
represented by 

<x i.j ' Vi ,j )/4e 


5 i, j 


and 


y 5|i.J ■ (y i,j ■ JH.jl'AE • 

The n-derivatives at these boundary points 
At j - JMIN, 


are represented as 


(3.47 a, b ' 
follows . 


x n [ i . j ‘ (x i,j»i x i,j )/4 ' 1 

y h|l,J " <y i,j.i - y i,j)/An , 
and at j - JMAX, 

X n|i,j " (x i,j ' x i,j-i )/An * 

y n|i,j “ (y i,j ' y i,j-i )/An . 


(3.47 c-f) 
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CHAPTER 4 


RESULTS AND DISCUSSION 


The analysis developed in the present study and described in the 
preceding two chapters is employed to solve the flow in constricted 
channels and several cascade configurations. The channel configurations 
considered here are the straight channel, which is employed primarily to 
verify the analysis developed, and a channel with an exponential 
constriction. The latter configuration is shown in Fig. 4a; its lower 
boundary is represented by the relation 


b, lower 


Ci exp [- ( 


(4.1) 


where the subscript b denotes the boundary, x^ is the x location where 

the maximum constriction is situated and Cj and C 2 are constants 
controlling the maximum height and the extent of the constriction, 
respectively. The equation for the upper wall of the channel is written 
as 

y =■ 1 - y . . (4.2) 

J b, upper ■'b, lower 


The cascade configurations considered in the present study are the 
flat-plate cascade and cascades with exponential, parabolic and 
Joukowski airfoils. These configurations are shown in Figs. 4b-e. The 
exponential airfoils are obtained using Eq. (4.1) for x^g £ x £ x^g, 

where x. „ and x_„ correspond to the x locations of the leading and 
Lh rh 

trailing edges, respectively. 

The parabolic arc airfoil is generated by the following equation: 
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X TE * X S X 1 


(4.3) 


y - C! + C 2 (x - ) , 


The Joukowski airfoil is generated by the following transformation. 



where 

z « x ♦ ly 

and 5 • C, + C s . (4.3 a-c) 

Here, C 3 and C„ are real constants and C s is a complex constant. 

The parameters C 3 and C 5 control the maximum thickness and the camber of 
the airfoil. As only symmetric airfoils are considered in the present 
study, C 5 is real. 

The results presented in this chapter have been obtained employing 
a H-grid or a channel-type of grid. A simple H-grid is the sheared 
Cartesian grid in which £ = £(x) and n = n(x,y). 


4.1 Resolution of Spatial Length Scales 

In viscous flows, the flow variables vary rapidly near the walls. 
To resolve these high gradients, a fine computational mesh is required 
near these boundaries. A non-uniform mesh i3 most suitable as it can 
provide a fine mesh in regions of high gradients without unduly 
increasing the total number of mesh points. 

The scalings obtained from the asymptotic analysis of Stewartson 
[38] indicate the order of resolution required for strong-interaction 
problems. Accordingly, streamwise mesh spacings should be of the order 
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of magnitude of Re and normal mesh spacings should be of the order 

-5/8 

of magnitude of Re . In the grids employed in the present 
calculations, at least five computational points are maintained within 
these length scales in regions of separation and near the trailing edges 
for cascade flows. The coordinate transformations employed to meet 
these mesh requirements are discussed next. 

In the streamwise direction, variable mesh spacing Ax^ is obtained 

using a geometric series for Ax^ and the resulting transformation can be 


written in a parametric form as follows: 


. (r 1 " 1 - 1) 

x i ' 4x * (r - 1) 


C = (i-1 ) AC , i = 1 , ... , ILE 


(4.5) 


Here, r is the stretching ratio, Ax, the mesh spacing at the inflow 
boundary, AC the uniform computational mesh spacing and ILE is the 
streamwise index corresponding to the leading edge location. Equation 
(4.5) is used with different values of r in various regions such as the 
airfoil surface and the wake, in order to obtain the required physical 
mesh spacings. 

In the n-direction, the following analytical transformation is 
employed: 

n = b + a tan . (4.6) 

c 

The parameters a and c in Eq. (4.6) provide control over the grid- 
point distribution. The constants a, b and c are obtained from the 
following conditions. 
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b, lower' 


n 


S, a prescribed value. 


At y - y 
At y - y 


n . and n 
min y 


n - n. 


(4.7) 


b, upper max 

The slope S controls the spacing near the boundaries. The transformation 
given by Eq. (4.6) is used to provide identical clustering near both 
boundaries, given by y - y p> lowep and y - y p> upper . Hence, b - 0.5. 

A typical grid is shown in Fig. 5. Here, the grid clustering in 
both the streamwise and normal directions has been reduced to improve 
clarity of the presentation of the coordinate lines in the figure. 


4.2 Results for Flow in a Straight Channel - Validation Study 

The analysis and the numerical procedure developed in the present 
research are first tested via a model problem of flow in a straight 
channel. For this purpose, a fully developed flow profile was 
prescribed at the inlet boundary and the pressure field was initialized 
so as to yield the known streamwise pressure gradient for a fully 
developed incompressible channel flow. This streamwise pressure 
gradient is given as 



The velocity and the temperature profiles at the inlet were 
obtained not from their known analytical expressions but by numerical 
solution of the governing equations for fully developed flow. This 
ensures that, if the algorithm is formulated in a consistent manner, the 
application of the procedure should recover the a fully developed flow 
throughout the entire channel. 
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This testing procedure was applied to the straight-channel 

configuration with metrics varying in the streamwise direction. From 

this study, it was found that the metric terms 5 and £ , which are 

x y 

associated with the streamwise pressure-gradient term p , should be 

S 9 * 

forward differenced. Also, it was found necessary that the streamwise 
pressure-gradient term should no longer be in the SCL form. The proper 
form of the pressure-gradient terms should be as given in Eq. (3.21) in 
Chapter 3- 

4.3 Results for Channels with Exponential Constriction 

The geometry of this channel is represented by Eq. (4.1) and is 
shown in Fig 4a. Results are obtained for three values of the ratio 
t/h of the maximum constriction to the channel width, viz., 0.1, 0.16 
and 0.2, for Re = 1500. The Reynolds number is based on the average 
velocity and channel height at the inlet station. As mentioned in 
Section 2.4.1, the inflow boundary conditions correspond to a fully 
developed flow in a straight channel. 

Figure 6a shows the distributions of the wall-pressure variable p^ 

and the wall shear parameter t w for the case of t/h =0.1. The wall- 

pressure variable p^ is defined as the difference between the pressure 

at the wall at a streamwise location i and that at the inlet boundary 
and can be written as 

p b ' p „, inlet ' p w, 1 • f 4 - 9 a > 

The wall shear parameter is defined as 
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wall 


(4.9 b) 


3U 


tgt 


w 

where U 


3n 


tgt 


is the velocity tangential to the body surface and n is the 


direction normal to the body. From Fig. 6a(i and ii), it is clear that 


P b varies linearly in the straight portion of the channel and behaves 


approximately similar to the body surface shape in the region of the 

constriction. From Fig. 6a(i), the pressure loss Ap across the 

constriction can be obtained to be 0.00875. From figures 6a(iii and 

iv), it can be seen that, far upstream and downstream of the 

constriction, t is asymptotic to the value of 6, corresponding to a 
w 

fully developed flow. The tendency of the flow to separate downstream 

of the constriction is indicated by the wall shear parameter approaching 

the value of zero. The grid employed for this case consists of 

(141 x 61) points. The value of Ax . is 0.008 which occurs at 
^ min 

locations upstream anc downstream of the maximum constriction where the 
flow is anticipated to separate. The minimum mesh step Ay mln employed 


in the y direction in the straight portion of the channel is 5.7 x 10 h 
Figures 6b and 6c show the wall pressure and the shear-parameter 
distribution for t/h of 0.16 and 0.2. From Fig. 6b(i), the pressure 
loss across the constriction can be obtained as 0.02, approximately. 

It is clear that the flow has separated downstream of the constriction, 
as indicated by the negative values of in Fig. 6b(iv). The results 

corresponding to the case of t/h ■ 0.2 are shown in Fig. 6c. Comparison 
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of Figs. 6b and 6c shows that the extent of separated region increases 
as t/h is increased from 0.16 to 0.2. 

The pressure loss across the constriction, some information about 
the grid and the location of the inflow and outflow boundaries with 
respect to the constriction are shown in Table 1 for the three cases 
discussed above. One of the important factors in obtaining the above 
results is the appropriate location of the inflow and outflow boundaries 
with respect to the maximum-constriction location. These locations are 
obtained by numerical experimentation, so as to ensure that the flow 
near these boundaries is nearly fully-developed. Accordingly, as shown 
in Table 1, for the case of t/h = 0.1, the inflow boundary should be 
located at least 5 channel heights upstream of the maximum constriction 
and the outflow boundary at least 1 8 channel heights downstream of the 
position of the maxiimum constriction. Moreover, these boundaries must 
be moved further away from the position of the maximum constriction as 
the ratio t/h increases. It is also observed that the pressure loss 
across the constriction increases with increase in t/h. 

Figure 7 shows the streamwise velocity profiles across the channel 
at various streamwise locations for the case with t/h « 0.2. From the 
enlarged-scale profiles shown in Fig. 7b, a small region of reverse flow 
can be observed. 

Through this study of the flow in constricted channels, the 
analysis has been tested for a variable cross-section channel and for 
separated flows. The regular behaviour of the flow solution in the 
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presence of separation establishes that the IPNs formulation includes 
the appropriate strong upstream interactions occurring in this flow. 

4.4 Results for Flat-Plate Cascades 

The configuration for this cascade geometry is shown in Fig. 4b. 

The primary purpose of examining these flows was to test the implicit 
implementation of the periodicity boundary condition for cascade flows 
and to test the validity of the IPNS formulation in the strong- 
interaction flows around sharp (cusped) leading and trailing edges. 
Results are obtained for various values of Re ranging from 1500 to 
16,000 and are shown in Fig. 8a-e. All of these results are obtained 
employing a grid containing (186 x 71) points. The point distribution 
is such that Ax is 0.01 and Ax_„ is approximately 0.005. The step 

L L 1 h 

size in the y-direction at the wall is 1.225 x 10 . As mentioned in 

Section 2.4.1, the inflow boundary conditions for this configuration, 
consist of the specified values of total pressure and total temperature, 
while static pressure is prescribed at the outflow boundary. 

Figure 8a(i) shows the distribution of the pressure p fe and the 

wake-centerline velocity. The corresponding behaviour of p w is inferred 

easily through Eq. (4.9 a). As the leading edge is approached, the 
pressure p^ rises to its maximum value and drops rapidly immediately 

downstream thereof. There is a gradual pressure drop in the flat-plate 
region, until the trailing edge, where a sharp drop occurs. Thereafter, 
the pressure rises smoothly to approach the prescribed value at the 
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outflow boundary. The wake-centerline velocity shows a smooth 
streamwise variation and is asymptotic to a value of unity near the 
outflow boundary. The wall shear, shown in Fig. 8a(iv), attains a 

maximum at the leading edge in an almost singular manner while 

exhibiting a small peak near the trailing edge of the cascade. The 
behaviour near the trailing edge is due to the sudden change in the 
boundary condition. Along the wall, the velocity component u was 
prescribed to be zero and along the wake centerline, u^ is zero. Hence, 

there is a discontinuity in u^., which is the primary contributor to the 

wall shear parameter. Through all this nonlinear behaviour, including 
that due to the sudden changes in boundary conditions at the leading and 
trailing edges, the solution is Quite regular, confirming again that 
upstream influence is appropriately included in the IPNS model. 
Qualitatively, similar behaviour is observed for Re = 3100 to 16,000, as 
shown in Figs. 8b-e. 

Figures 9a and b show the distribution of p^ and x^ for various 
Reynolds numbers. With increase in Re, the peak in p^ at the trailing 
edge decreases, whereas that in t w increases. 

14.5 Results for Cascades of Exponential Airfoils 

The exponential airfoil cascade geometry is as shown in Fig. 4c. 
The computational grid, the point distribution and, hence, the mesh 
sizes employed are the same as for the flat-plate cascade configuration 
discussed in the previous section. 
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4.5.1 Effect of Thickness 


Figure lOa-c shows the distributions of p and t for an 

D W 

exponential airfoil cascade for Re - 1500 and thickness ratios (t/c) of 
0.05, 0.075 and 0.1. Near the leading and trailing edges, the behaviour 
of pressure is similar to that for the cascade of finite flat plates. 

In the region of the airfoils, as the overall flow accelerates up to the 
maximum thickness location, the pressure falls to a minimum, while the 
shear parameter increases to a local maximum. The pressure and wall- 
shear parameter distributions for various values of t/c are shown in 
Figs. 11a and b, respectively. Superimposed on these figures are the 
corresponding flat-plate cascade results, i.e., the case of t/c =0. 

From Fig. 11a, it is clear that the maximum p^ increases with increase 

in t/c, while the peak at the trailing edge diminishes. The value of t 

w 

at the maximum-thickness location also increases with increase in t/c. 
The tendency for the flow to separate downstream of the maximum- 
thickness increases with increase in t/c. This can also be inferred 
from the pressure gradient, which becomes increasingly adverse in this 
region as t/c is increased. Also, the extent and intensity of the 
separated region are reduced, in comparison with the corresponding 
channel flow. This is because of the upstream influence of the higher 
velocity of the fluid downstream of the TE as compared to that of the 
fluid downstream of the constriction in the channel with zero slip at 
the walls. 
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4.5.2 Effect of Reynolds Number 


Results for various values of Re between 3100 and 15,000 are shown 
in Fig. 12a-d. The qualitative behaviour of the flow properties 
examined for this range of Re is similar to that for Re = 1500, 
discussed in the previous section. Figure 13 shows the effect of Re on 
the pressure and wall-shear distribution. From Fig. 13a, it can be seen 
that the magnitude of p fe at the maximum-thickness location decreases 

with increase in Re. The pressure peak near the trailing edge also 
diminishes with increase in Re. It can also be seen that the presoure 
gradient downstream of the maximum thickness becomes increasingly 
adverse with increase in Re. A similar behavior in pressure can be 
seen, to a smaller degree, upstream of the maximum thickness, xhe 
possibility of flow separation in these adverse pressure-gradient 
regions is evident from the distribution of t w shown in Fig. 13b. In 

fact, for the case of Re = 15,000, a small separated region is present, 
as indicated by the negative values of t^. It is also evident unat the 

separation region is centered around a point situated downstream of the 
maximum-thickess position and upstream of the trailing edge, i.e., as Re 
is increased, the separated region extends in both directions about this 
point. The analyses of H.K. Cheng and F.T. Smith [39] and Smith, 
Stewartson and Kaups [40] show similar results around a cusped trailing 
edge. 
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4.5.3 Effect of Grid Refinement 

The effect of grid refinement is important in order to establish 
the reliability of the solutions obtained, especially when there are no 
experimental results with which to compare. Therefore, in the present 
study, the streamwise grid is refined in the region downstream of the 
maximum-thickness location, in order to study its effect on the 
separated region. For all the cases discussed in the previous section, 
the streamwise step Ax fc (i.e., Ax near the maximum-thickness location) 

was approximately 0.045 and the streamwise step size decreased gradually 
to a value of 0.005 near the trailing edge. In the grid-refinement 
study, a finer streamwise grid in this region was obtained by ensuring 
that AXj. = 0.025 and Ax^. r = 0.005, and employing a computational grid 

consisting of (191 x 71) points. 

A comparison of the results obtained employing the two grids, for 
Re - 1500 and 15,000, are shown in Figs. 14a and 14b, respectively. The 
wall-shear distribution agrees well and, in the case of Re = 15,000, the 
extent of the separated region remains almost the same, as seen from 
Fig. 1 4b < i i ) . Grid refinement has a slightly more significant effect on 
the pressure distribution, particularly for the case of Re = 15,000. In 
that case, the peak near the trailing edge diminished with grid 
refinement. 

4.6 Results for Cascades of Parabbl ic-Arc Airfoils 

A schematic of this cascade configuration is shown in Fig. 4d. All 
the results discussed in this section are obtained employing a grid 
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consisting of (186 x 71) points, with the streamwise mesh sizes near the 
leading and trailing edges being 0.01 and 0.005, respectively. 

H.6.1 Effect of Thickness 

Results in the form of the distribution of p fe and t w are shown in 

Figs. 1 5a-c, for t/c of 0.05, 0.075 and 0.1. Figure 15a(ii) shows the 
pressure distribution and the development of the wake- centerline 
velocity for the case of t/c - 0.05. An important difference between 
the pressure distribution for an airfoil with a wedge-shaped trailing 
edge such as the parabolic-arc airfoil and that for an airfoil with a 
cusped trailing edge such as the exponential airfoil, is that at the 
trailing edge the peak in the p fe di stribution for the exponential airfoil 

cascade is replaced by a slope discontinuity for the parabolic-arc 

airfoil cascade. Also, from Fig. 15a(iv), it can be seen that the peak 

in the x distribution at the trailing edge has diminished compared to 
w 

that at a cusped trailing edge. From Fig. 16a, it is seen that p fe at the 
maximum-thickness location increases with increase in thickness ratio 

and the pressure gradient on the downstream side of maximum-thickness 
location becomes increasingly adverse. The wall-shear distribution 
exhibits a cross-over, with a small separated region near the trailing 
edge for the case of t/c ■ 0.1. 


63 



4.6 .2 Effect of Reynolds Number 

Figures 17a-d show the results for various Re ranging from 3100 to 
15,000. From the velocity distribution in Fig. 17c(ii) and the wall- 
shear distribution in Fig. 17c(iv), corresponding to Re - 11000, it can 
be seen that a separated region exists, starting at a point upstream of 
the trailing edge, with the flow reattaching at a point in the wake. 
Further increase in Re results in the movement of both the separation 
and reattachment points away from the trailing edge. Therefore, for 
flows about wedge-shaped trailing edges, the separated region is 
centered at the trailing edge. This is due to the continued 
deceleration of the fluid up to the wedged TE whereas, for cusped 
trailing edges, the deceleration decreases as the fluid approaches the 
cusped TE. This result is also supported by the analysis of Cheng and 
Smith [39]. Comparison of the pressure and wall-shear distributions for 
the various values of Re discussed above is shown in Fig. 18. 

4.6.3 Effect of Mach Number 

In all of the results discussed so far, was approximately 0.008. 

Results have also been obtained for Mach number ranging from 0.035 to 
0.49 and are shown in Figs. 19a-c. All of these results are obtained 
for Re ■ 15,000. Figure 20a shows that the effect of compressibility, 
resulting from the increase in M^, becomes apparent on the wall-pressure 

distribution only for the case with - 0.49. But the wall-shear 

distribution as well as the extent of separation remain almost unchanged 
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for all Mach number values considered. This result is due to the fact 
that the strong trailing-edge singularity overwhelms the effects due to 
the Mach number change, in the range considered. Figure 21 shows the 
contours of the static pressure p, for M = 0.49 and Re = 15,000. 

^•6.4 Effect of Grid Refinement 

In order to ensure the accuracy of the behaviour of the flow near 
the wedge-shaped trailing edges discussed above, it was considered 
necessary to refine the grid in this region. Hence, results were 
obtained employing a grid consisting of 231 points in the streamwise 
direction, so that the Ax,^ = 0.001 . Results obtained using this grid 

are compared with those obtained with the grid of (186 x 71) points as 

discussed above, for Re = 3100 and 11,000, and are shown in Figs. 22a 

and 22b, respectively. The shear parameter t remains unaffected by 

w J 

this grid refinement for the two values of Re considered hers. Although 
there is a slight shift in the level of the pressure as a result of 
refining the grid, the behaviour near the trailing edge remains 
unchanged. 

^ ^ Results for Cascades of Joukowski Airfoils (Modified Leading Edge) 
The cascade conf igurations considered in the previous sections were 
made up of airfoils with either cusped or wedge-shaped leading edges. 

The procedure developed in the present study was applied next to cascade 
passages with rounded leading edge blades such as the Joukowski airfoil, 
for which a typical configuration is shown in Fig. 1. This 
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configuration presented some difficulties to the numerical solution in 
the leading-edge region. These difficulties were traced back to the use 
of the sheared Cartesian grid for these cascades. It was observed that 
this type of grid exhibits a large discontinuity in the coordinate- 
transformation metrics at the leading edge. Hence, the leading-edge 
region was modified by replacing it with a wedge, as shown in Fig. 4e. 
Results have been obtained successfully for the modified configuration, 
even with the sheared Cartesian grid, for Re » 150 and 310 , and are 
shown in Figs. 23a and 23b. The computational grid employed consists of 
(121 x 61) points and the value of the ratio t/c for the airfoil is 

0 . 068 . 

Figures 23a(i) and 23b(i) show the p b distribution and the 

development of streamwise velocity along the wake, for Re =150 and 310, 
respectively. It is evident that, as Re is increased, the peak in the 
P b distribution near the trailing edge diminishes, as would be expected 

for flows over cusped trailing edges and as observed for the 
exponential-airfoil cascade. The wall-shear distribution, shown in 
Figs. 23a(ii) and 2 3b ( i i ) , has a slope discontinuity near the maximum- 
thickness location. This is primarily a reflection of the discontinuity 
of the slope of the body surface where it changes from the wedge to the 
Joukowski-airfoil profile. 

The difficulty in resolving the flow near the rounded trailing 
edges with an H-type grid arises mainly due to the orientation of the 
grid in this region. To circumvent this problem, a body-oriented grid 
such as a C-grid is more appropriate for rounded leading-edge blades and 
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should be employed. Therefore, a major effort was directed, in the 
present study, to develop a hybrid C-H grid which judiciously employs 
both C- and H-grids, thereby benefiting from the advantages of both of 
these types of grids. The generation and implementation of the C-H grid 
for a cascade of airfoils with rounded leading edges is discussed in the 

next two chapters. 

4.8 Convergence Study 

A typical convergence history of the solution is shown in Fig. 24 
in terms of the maximum absolute error and the relative r.m.s. error in 
pressure, for the case of the parabolic-arc airfoil of thickness ratio 
0.05 and Re = 6300. From this figure, it can be seen that convergence is 
achieved in approximately 1 60 iterations. All the other cases discussed 
in the previous sections also required a similar number of iterations. 
The definition of the errors and the tolerances used in the present 
study are described in Section 3.^. In all the results presented thus 
far, 3 local iterations were performed at each streamwise location in 
order to update the non-linear terms. It was found that quasi- 
linearization about a previous streamwise location was insufficient, 
especially in regions of large streamwise pressure gradients, such as 
the leading- and trailing-edge regions. Finer streamwise steps alone 
cannot resolve these high gradients due to their near-singular 
behaviour. 

In order to study the effect of the time step At associated with 
the pressure gradient terra (Eq. 3.25), a simple case of the flat-plate 
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cascade was considered. A constant step size of Ax = 0.05 was employed 
for this case. The decay of the maximum absolute error and the relative 
r.m.s. error in pressure was studied for various values of At in the 
range between 0.03 and 0.1. The corresponding convergence histories are 
shown in Fig. 25. For all values of At used, the convergence behaviour 
is identical for approximately the first 20 iterations. Beyond this 
stage, the convergence rate increases as At is increased from 0.03 to 
0.06 and decreases with further increase to a value of 0.1. For the 
value of At = 0.06, a total error-reduction of approximately 6 decades 
is achieved in 200 iterations. Hence, it is observed that the procedure 
is rather sensitive to the value of the time step At used. 

The convergence behaviour for a flat-plate cascade with variable 
streamwise mesh size was also studied. For this case, the mesh width 
was varied according to Eq. (4.5). Also, At was varied in the 
streamwise direction, according to the following relation: 

At i = C (x i+1 " x i } * (4.10) 

Here, C is a constant and subscript i denotes the streamwise location. 

Figure 26 shows the convergence history for various values of the 
parameter C. It is clear that for the first 40 iterations, the 
parameter C does not affect the convergence rate. Thereafter, an 
increase in C results in faster convergence rates, for the range of C 
considered in the present study. 

Finally, a convergence study was performed for the cascade with 
exponential airfoils with t/c = 0.05 and Re = 1500. The grid employed 
for this study consisted of (191 x 71) points. Figure 27 shows the 
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convergence history for values of C ranging from 2 to 5. As C is 
increased from 2 to 4, the convergence rate improves, but further 
increase in C results in a slower convergence rate. It can be seen that 
the slope of error curve for C = 3 is steeper than that for 4 so that 
C = 3 may be considered as a near-optimum value. Then, with AtVAx.. = 

3, results were also obtained for Re = 25,000 for a cascade of 
exponential airfoils and are shown in Fig. 28. 
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CHAPTER 5 


GENERATION OF HYBRID C-H GRID 

5.1 Introduction 

The success of a numerical method for analyzing viscous flows in 
complex configurations depends on the proper choice of a coordinate 
system. It is usually desired that these coordinates be aligned with 
the problem boundaries. Also, as the governing equations of motion in 
the present study (see Chapter 2) involve a parabolizing approximation 
in the streamwise direction, it is important that the coordinate system 
employed be body-oriented and near-orthogonal. 

For the study of subsonic viscous flows around isolated airfoils, 
two types of grids are often used. These are the O-grids and the C- 
grids. For analysis of flow through cascade passages formed by 
airfoils, the H-grid or the channel type of grid and the C-grid 
configurations are commonly used. Sockol [41] and Dulikravich [42] have 
generated C-type of grids for cascades. A typical C-grid is shown in 
Fig. 29a. Steger et al . [8] have also employed C-grids for calculations 
of viscous and inviscid flows through turbomachinery cascades. The C- 
type of grids are especially suited for subsonic viscous flow 
calculations, as they provide good resolution near the rounded leading- 
edge regions of the blades. From the results described in the previous 
chapter for flow through a cascade of Joukowski airfoils, it is clear 
that the channel type or H-grid is inappropriate, particularly near the 
rounded leading edges. A C-type of grid is satisfactory in most regions 
of the cascade flows, except further upstream of the stagnation point on 
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the blades. In this region, the grid density decreases rapidly with 
increasing distance upstream of the stagnation point. A possible means 
of minimizing this difficulty is to employ a channel or H-grid in the 
upstream region. Eiseman [ 4 3 ] has employed Cartesian extensions to an 
O-type of grid for cascades, in the regions upstream of the leading 
edges as well as downstream of the trailing edges. 

Channel grids have been widely used for cascade-flow analyses. 

Chima and Johnson [10] have employed the H-type of grid for solving the 
Euler and NS equations through a cascade of bicircular-arc airfoils. 
Channel grids are easy to generate and can be conveniently aligned with 
the inflow. Also, implementation of the periodicity boundary condition 
is easy for cascades with low stagger angle. In this case, a grid line 
emanating from a particular point on the lower boundary ends at the 
corresponding periodic point on the upper boundary. With this grid 
distribution it is also possible to impose the periodicity con.-* -ion 
an implicit manner (Ref. [44]). This H-grid distribution is termed 
-line-periodic' and is shown in Fig. 29b. For cascade configurations 
with large stagger or when the blades have rounded leading edges, the 
line-periodic H-grid becomes highly skewed and non-orthogonal . If this 
H-grid is modified to minimize the non-orthogonality, the implicit 
imposition of periodicity has to be sacrificed. A typical modified H- 
grid is shown in Fig. 29c and is called a 'region-periodic' grid. In 
this type of grid, the point distribution along the upper and lower 
boundaries is periodic, but corresponding periodic points are not 
connected by the same coordinate line. 
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Therefore, a suitable grid for cascade nous is one that combines 
the advantages of both the C- and H-grids while minimizing their 
disadvantages. This results in a multi-rectangular computational domain 
in terms of the transformed coordinates. The multi-block structure and 

the generation of such a grid are described in the following sections of 
this chapter. 

5-2 Multi-Block Structured Grids 

The flow region in the physical domain can be subdivided into 
several sub-regions, depending on either the geometrical complexity of 
the configuration or the need for different types of grids in various 
regions, as mentioned in the previous section. These sub-regions are 
termed ’blocks’. Lee et al. [45] have generated a 3 -D body-fitted 
coordinate system around a wing-body-nacelle configuration by dividing 
the computational domain into multiple rectangular blocks. A single 
rectangular block in the computational domain would be insufficient to 
resolve the various physical corners of this complex configuration. 
Multi-block structured grids have also been employed for the complex 
geometries of the tri-element augmentor wing by Sorenson [46] and the 
breaking surface wave by U. Ghia et al. [47]. Also, Coleman [48] has 
employed multiple segments for generating alternate grids for flow over 
single airfoils. The 0-type of grid with Cartesian patches for a 

cascade generated by Eiseman [4 3 ] can also be viewed as a multi-block 
structured grid. 
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Multi-block structured grids can be broadly classified into two 
types. The first type consists of 'patched' or 'zonal' grids where the 
various blocks share a common boundary called the interface. The other 
type of grids consist of 'overlaid' or 'overset* grids. In this case, 
the various blocks do not have a common boundary. The idea of overset 
grids has been pursued by Steger [^9] in his development of the Chimera 
grid, which has also been extended to three dimensions by Benek et al. 
[50]. A typical overset grid is shown in Fig. 30a. The solution of the 
flow equations employing this type of grid requires additional transfer 
of information across the boundaries of various grids. 

Patched grids can be further classified into joint and disjoint 
types of grids. The disjoint patched grid is one where the family of 
the coordinate lines crossing the patch interface are discontinuous 
across the interface, as shown in Fig. 30b. In the joint type of grid, 
shown in Fig. 30c, the coordinate lines are continuous across the 
interface, although their slope may be discontinuous across the 
interface. These gradient discontinuities are present when the 
coordinates in each of the sub-regions are generated separately and 
patched together. Some of these discontinuities can be relieved by 
suitable modification of the boundary-point distribution along the 
interface, followed by regeneration of the coordinates in the regions 
affected by this redistribution. These slope discontinuities can be 
eliminated by generating the coordinates in all the sub-regions 
simultaneously, without specifying the boundary-point distributions 
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along the interface. The disadvantage associated with a such a 
procedure is that the computer program can become quite complicated, 
depending on the number and orientation of the blocks employed to define 
the complex physical configuration. The latter direct multi-block grid 
generation concept has been employed by U. Ghia et al. [21] and Coleman 
[48]. The composite grid-generation procedure developed by U. Ghia et 
al. [ 21 ] is employed in the present study and the details are discussed 
in the next two sections. 

5.3 Computational Domain for Hybrid Cascade Grids 

For a general staggered cascade, the physical domain can be divided 

into sub-regions as shown in Fig. 31. The region BCD' B' is the portion 

of the C-grid adjacent to the upper surface of the lower blade. The 

regions BIHG and E'D'EM are the channel or H-regions. The computational 

regions corresponding to these sub-regions are shown in Fig. 32. Points 

E" and B" are 'special-points' at which 5 sub-regions come together, 
c c 

'Special-points' commonly occur on the domain boundary or on the 
interfaces between sub-regions in a multi-block grid system. These 
points either have a non-standard number of immediate neighbours when 
they are the vertices of a computational cell or have a non-standard 
number of faces when they are the cell-centers, as illustrated in Fig. 
33. 

For the ease of visualization, the multi-block structured grid in 
the 2-D physical plane is best represented by a 3”D surface in the 
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computational domain. This 3-D surface is obtained by bringing the sub- 
regions in Fig. 32 together, along the arrows. The subregions BCO’B" 

and IBB-L are bent out of the plane of the paper at the special-point 
c 

. „ uvrtir r.p.wp ,, B and C B^E^E* are bent at 
B". Similarly, the sub-regions HKG G, GG ana c 


E ». The resulting computational domain is shown in Fig. 34. For a 
cascade with zero stagger, the sub-regions G E c B B c and 

B-b'D-G’ are not needed. The sub-region boundaries in the physical 
c 

domain for the unstaggered grid are shown in Fig. 35. The solution 
procedure for generating the coordinate system is discussed in the next 

section. 


5.14 Solution Procedure 

The numerical transformation procedure consists of determining the 
boundary-oriented coordinates as the solutions of the following Poisson 
equations : 

(5.1 a) 

V 2 S = P 


and 

(5.1 b) 

V 2 n = Q 

where V* is the Laplacian operator in Cartesian ooordinates, 

( 5 ,n) are the transformed surface-oriented coordinates and 
P and Q are the forcing functions used to provide control on the 

coordinate clustering and orthogonality. 

The boundary conditions used for the governing equations (5.1) are 
of the Oirichlet type. These consist of the prescribed values of x and 
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y along the boundaries r = r an rt c anH _ 

^in «max and " ' U and \ax* In 

order to Implement these boundary conditions, the roles of the 

independent (x.y) and dependent <5,n) variables have to be interchanged. 

Detailed analysis of this inversion is given by U. Ohia and K. Ghla [51] 

and the resulting 'inverted' equations can be written as follows: 

3 X « * 2b x 5n * 0 x nn * j2 <P * 0 *„> - 0 

and 


a * 2b y Cn + C y nn + J 2 (P y^ + Q % ) = 0 

where a = x 2 + v 2 

n J n ’ 


b ‘ ' <X e x n * h V 


° ' X 5 * h 


and J=x„y — x y 
t n n 5 


(5.2) 


In the present study, the alternating-direction implicit (ADI) 
method is employed to solve the governing equations ( 5 . 2 ). The 
computational domain for the and the ^-implicit sweeps of this 
procedure can be obtained from Fig. 31 and are shown in Figs. 36a and 
36b, respectively. In Fig. 36 , points denoted by alphabets with 
superscript * or -, refer to the neighbouring points. It should be 

noted that the subregions IL~Bj B - , hkVg" and cVe-B^ are encountered 

twice in the n-implicit sweep and, hence, the solution in these regions 
at the end of this sweep is already at the <n*1 ) level of the ADI 
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procedure. Therefore, these sub-regions must be excluded from the 

computational domain in the ^-implicit sweep. 

The treatment of various derivative terms that appear in Eqs. (5.2) 
needs special considerations around the points B" and E”. In the 

present study, these special points are considered to be the center of a 
non-standard, 5 -sided cell. The selection of the appropriate 
neighbouring points involved In the discretized form of the governing 
equations at the vertices of the 5 -sided cells is detailed in 

Appendix C. 

5.5 Typical Grids 

Figure 37a shows the hybrid C-H grid generated for a staggered 
cascade using the composite procedure just described. This grid is 
region-periodic and consists of 1620 points. The coordinates are 
uniformly spaced at the inflow boundary and are clustered along the 
wake. Figure 37b shows a hybrid grid for an unstaggered cascade of 
Joukowski airfoils. This grid is line-periodic and consists of 7103 
points and is employed for obtaining a flow solution also using a 
composite solution procedure. This procedure for calculating the flow 
through a cascade employing a hybrid C-H grid is discussed in the 

chapter. 
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CHAPTER 6 


determination of f low through JOUKOWSKI CASH APP- 
OSING HYBRID C-H RRTn 

In the preceding two chapters, the inadequacy of the H-grid or the 
channel type of grid to resolve the flow around rounded leading edges, 
was pointed out. Also, the difficulties of using C-grids for cascade 
flows were discussed simultaneously. The disadvantages of using these 
grids in certain regions of cascade flows were observed. Hence, a 
viable alternative grid is the hybrid C-H grid. This hybrid grid, with 
its associated multiple rectangular computational domain, was developed 
in the preceding chapter, for use with cascade-flow studies. Analysis 
of flow over complex geometries, using multi-block structured grids, has 
been performed by many researchers. Benek, Steger and Dougherty [52] 
have employed overset grids to obtain the solution of the Euler 
equations about a supercritical flapped airfoil. They observed that the 
solution in the transonic regime exhibited an ill-defined shock at the 
grid boundaries. Eberhardt and Baganoff [53] developed characterestic 
boundary conditions, which alleviated the above mentioned difficulty. 
Norton, Thompkins and Haimes [18] have employed a patched grid of the 
joint type, consisting of sheared grids and O-grids, and a cell-centered 
implicit scheme to solve the complete Navier-Stokes equations for flow 
in turbine cascades. Rai [19] has used a system of patched and overlaid 
grids and obtained the solution of the thin-layer Navier-Stokes 
equations for a rotor-stator combination. When such a grid is employed, 
the solution on one grid system at the patch boundary has to be 
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interpolated and the boundary condition for the other grid system should 
be obtained in a manner that conserves all fluxes crossing the interface 
boundaries. In the present study, the semi-elliptic analysis developed 
in Chapter 2 is employed to determine the flow through a cascade of 
Joukowski airfoils, using the hybrid C-H grid described in the preceding 
chapter. The solution procedure for the governing equations of motion 
on the hybrid grid is described next. 

6 . 1 Computational Procedure 

The solution procedure used with the hybrid grid is, in principle, 
similar to that for the H-grid and consists of two steps, as described 
in Section 3.2. For a hybrid C-H grid, these two steps are represented 
schematically in Fig. 38. The regions denoted by HI and H2 correspond 
to the H-grid regions of the hybrid grid and the regions Cl and C2 
correspond to the two C-regions. The first step of the solution 

procedure, during which the flow variables Q are updated from time level 
n to (n + 1/2), is denoted by a thin arrow; the second step, in which 
the pressure field is updated from time level (n+1/2) to (n+1), is 
denoted by a thick arrow. The computation proceeds from the inflow 

boundary of region HI, with Q being updated first in this region. The 
boundary conditions used at the upper and lower boundaries in this 
region are the periodicity conditions, described in Section 2.4.3. 

Next, the first step of the solution procedure is applied to region Cl. 
In order to employ a marching type of procedure, conditions should be 
prescribed along a line such as AB, in Fig. 38b. The flow conditions 
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along AB are obtained from the solution along the corresponding periodic 
line E 1 F 1 at the time level n. When the first step of the solution 

procedure is applied to region C2, the flow conditions along could 

be obtained from the latest available solution along the line EF. But, 
in the present study, in order to retain symmetry, even in the transient 
stages of the solution procedure, the flow conditions along a B are 

obtained from the solution along EF at time level n. The boundary 
conditions in the C regions are the no-slip condition and the prescribed 
temperature at the airfoil surfaces, while boundary conditions at the 
interface are obtained from the solution in region HI. These interface 
boundary conditions are discussed in a later section of this chapter. 
Next, the solution in region H2 is obtained by marching, knowing the 
solution along the outflow boundaries of region HI and the two C 
regions. 

The second step of the solution procedure, during which the 
pressure field is updated, proceeds from the outflow boundary. The 
outflow boundary condition on pressure is applied directly during this 
step. Using Eq. (3.3”), the pressure field is updated from a time level 
(n+1/2) to (n+1 ) in region H2 and then in the two C regions. This step 

of the procedure is then continued into region HI, knowing the values of 
n +1 

p along the boundaries of the C region at the C-H1 interface and the 
first computational line in the region H2. The flow variables at the 
inflow boundary are updated as discussed in Section 3 . 3 . The two steps 
are repeated until convergence is achieved. 


80 


6.2 Verification via a Mathematical Model Problem 

The correctness of the formulation and programming of the solution 
procedure on the hybrid C-H grid was checked through application to a 
model problem. The verification of the scheme by recovering the known 
solution for fully developed flow in a straight channel cannot be 
employed here, as the C-H hybrid grid is not appropriate for this 
configuration. Hence, a mathematical model problem is employed. In 
order to check the first step of the solution procedure, the system of 
governing equations, given by Eq. (3*23) is considered and can be 
written as 


* * 

M . |E . b 

35 3n 


3p 

3n 




( 6.1 ) 


The quasi-linearized form of this equation is obtained from Eq. 


(3.24) and is 


A A Q + A5 


{(B * 1 -B 1 ) A*Q } + b i + 1 f- U ^) 1 A*q] 
on v <3n 3 q 


*i -*i -i 

E - A Q 


A5 {I- (B* 1 Q 1 ) + b 3 
3n 


( 3P)i 

v 3n 


i+1 .i + 1 
- a A p 


Next, Q is taken to be an analytical function as follows. 

In regions HI and H2, 

Q = Q x sin (oil 5) cos (2ir n) + Q 2 » 
and, in regions Cl and C2, 

Q = Q 3 cos (a 2 n + I 2 ) ^*3 a * b ) 
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where 


Q*, Q 2 » Q* i a and 0 are specified vectors. 


In general, Q chosen as in Eq. ( 6 . 3 ) will not satisfy the governing 
equations (6.1). Therfore, the governing equation (6.1) is modified by 


including in it an extra source term, R(5,n), such that Q is a solution 

of the modified equation. This modified equation is written as 
* * 


9E + 3F 

H 3n 


b lE _ i_ (p ) 

9n 3n v' 


•ft. 


R(C.n) 


(6.4) 


The term R(5,n) could be evaluated analytically using Eqs. (6.3) 
and (6.1). In this case, the solution that would be obtained by solving 
Eq. (6.4) numerically will not exactly equal the chosen analytical 


function 0(5, n) and, hence, the correctness of the formulation may not 
be established accurately. It is, therefore, important that the source 
term R(5,n) be evaluated using the discretized form given by Eq. (6.2). 


This model problem was studied for various c^, a 2 , fj, and S 2 . The 
model problem of uniform flow when all the boundary and initial 
conditions correspond to uniform flow is a degenerate case of the model 


problem discussed above and corresponds to 5 - constant. In all the 
cases considered, the exact solution was reproduced. It was found that 

when Q was allowed to vary rapidly in the 5 direction, large errors were 
observed between the prescribed and computed values of the solution. 

This was due to the rather large truncation errors involved due to large 
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gradients in Q in the C direction. This implies that finer grids in the 
£ direction have to be employed to resolve these high gradients. 

6.3 Interface Boundary Conditions 

The implementation of the semi-elliptic formulation on a hybrid 
grid exposes two small portions of the interfaces of the grid. Hence, 
the boundary conditions on these interfaces have to be considered. 

The boundary conditions in the C-regions consist of the no— slip 
condition and the specified temperature at the walls, together with the 
conditions specified along the interface of the C and HI regions. These 
interface conditions are obtained from the flow solution in the HI 
region. As the governing equations are of order seven with respect to 
the n direction, M conditions are yet to be specified. These conditions 

T 

could be obtained from the solution, Q = (p, pu, pv, pe t ) , computed in 

region HI . But the physics of the flow dictates how disturbances travel 
in the flow and this must be incorporated in the numerical scheme used. 

A characteristics analysis for the quasi— one-dimensional Euler equations 
suggests that for subsonic flows only three conditions can be specified 
along the upstream boundary and that one condition is to be specified at 
the downstream boundary. This idea has also been employed by Eberhardt 
et al. [53] and Bush [20]. The former reference deals with the solution 
of the Euler equations while, in the latter, the complete Navier-Stokes 
equations are solved. In regions of flow where viscous effects 
dominate, the characteristic conditions obtained from the Euler 
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equations are not appropriate. However, these conditions should provide 
a reasonable approximation. In the present study, the three conditions 
used at the interface consist of specified values of the density p, the 
x-component of the velocity u, and the static pressure p. The 
additional condition needed to close the system is the reduced normal 
momentum equation (2.22) written at the half-mesh point off the wall. 

6.4 Application to Joukowski Airfoil Cascade 

The procedure described in the previous sections is applied to 
determine flow at Re = 150 through a cascade passage formed by Joukowski 
airfoils. The hybrid grid that was employed is shown in Fig. 37b. The 
extent of the C-region, measured in terms of the physical length of BC 
in Fig. 35, is approximately 1 0% of the airfoil chord. This value was 
arrived at from the results obtained previously employing an H-grid and 
corresponds to the nearest streamwise location upstream of the leading 
edge of the cascade where p^ and u are approximately zero. The 

streamwise distributions of the surface pressure p fa , the wake-centerline 

velocity and the wall-shear parameter t are shown in Figs. 39a and 39b. 

W 

The surface-pressure distribution is similar to that obtained for an 

exponential airfoil cascade, except that the behaviour near the leading 

edge is confined to approximately 1/4 of the chord. The t distribution 

w 

near the trailing edge exhibits a rise, as should be expected for low Re 
flows around cusped trailing edges. The velocity distribution along the 
stagnation line exhibits an overshoot in the C-region. This is probably 
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due to the fact that that the stagnation line represents a singularity 


in terms of Q and needs special treatment, paticularly for the y- 
momentum equation. Figure 39d shows the pressure contours superimposed 
on the hybrid grid. It is clear from this figure that there is an 
oscillation present along the first line of H2 region and this anomalous 
behaviour starts at the 5-sided cell. This oscillatory nature is due to 
the discontinuity in the y-component of velocity v, between C and HI 
regions. This discontinuity leads to a sudden jump in v^ which appears 

in the viscous terms. This is supported by the fact that with increase 
in Re, the oscillatory behavior in p vanished. The results for Re = 600 
are shown in Fig. 40. The pressure contours shown in Fig. 40c do not 
exhibit any oscillation, but the anomaly near the special cell persists. 
Figure 40e shows the development of velocity profile for Re - 600. Th- 
velocity profiles are well behaved, in spite of the anomaly in pressure, 
because these anomalies are of the order 0.01 . Figure 41 shows the 
convergence history for this Reynolds number. The e abs and e rel » 

defined in section 3-3, are reduced by approximately 4 orders of 
magnitude in 100 iterations. The anomalous behaviour in pressure is due 
to the discontinuity in the metric derivatives along the coordinate line 
containing the downstream face of the five-sided cell. Along this line, 
the chain rule conservation law (CRCL) form, instead of the SCL form, is 
employed. The results are shown in Fig. 42. The pressure contours 
shown in Fig. 42c no longer exhibit the unrealistic behaviour near the 
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five-sided cell. Effort is being made to resolve this problem near the 
stagnation line in order to obtain solutions at higher Re. 
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CHAPTER 7 


CONCLUSION 

A semi-elliptic formulation has been developed for the analysis of 
subsonic viscous flows. The governing equations were obtained from the 
NS equations by neglecting streamwise viscous diffusion terms, but 
retaining upstream interaction via appropriate treatment of the 

streamwise pressure^gradient term. 

A numerical procedure has been developed to solve the semi-elliptic 
equations. The discretization of the governing equations was performed 
such that the terms representing the inviscid contributions form a 
consistent set of equations by themselves. Appropriate discretized 
forms of the metric terms associated with the SCL form of the equations 
were obtained by requiring that a uniform flow solution be recovered in 
the interior when all boundary and initial conditions correspond to 
uniform flow. The discretized representation of the streamwise 
pressure-gradient term, as well as the metrics associated with it, was 
obtained by application of the procedure to the fully developed flow in 
a straight channel. Through this study, it was found necessary that the 
streamwise pressure-gradient term not be in SCL form and that the 
metrics associated with this term be discretized using forward 
differences. 

The procedure developed was applied to several flow configurations 
such as channels with exponential constrictions and cascades of airfoils 
of various shapes. The technique was demonstrated to be adequate for 
strong-interaction flows, where boundary-layer separation and/or sudden 
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changes in boundary conditions due to the presence of sharp 
leading/ trailing edges are present. A parametric study was carried out 
by varying Re for flows through cascades of exponential and parabolic- 
arc airfoils. This revealed a considerable difference between the flow 
behaviour near cusped and wedge-shaped trailing edges. The effect of 
Mach number was studied for a cascade of parabolic-arc airfoils in the 
range of between 0.008 and 0.49, and was found to be minimal on the 

shear parameter, but to be evident in the pressure distribution for the 
highest considered. The effect of grid refinement was also studied 

in order to establish the accuracy of the results obtained. In all of 
the calculations, the grids employed were reasonably fine and the step 
sizes near the TE and near the point of separation were chosen to be 
within the scalings of the triple-deck theory. The grids employed were 
considered to be adequate because the minimum step sizes were well 
within the triple-deck scalings. The grid independence of the results 
was also supported by the fact that the results remained unchanged with 
further refinement of the grid. As regards the convergence rate, for 
most of the flow solutions considered here, the maximum error was 
reduced by four to five orders of magnitude, within 150 iterations. 
However, the procedure remained sensitive to the magnitude and the 
spatial variation of the time step At. 

The application of the procedure to cascades with rounded leading 
edge regions revealed some difficulties. These difficulties are 
attributed to the use of an H-grid for such configurations. To resolve 
these difficulties, a hybrid C-H grid, which is more appropriate for 
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C - 3- 



The 


such cascades, was generated using a composite solution procedure, 
basic solution procedure was then extended for use with the multi-block 
structured computational domain corresponding to the C-H hybrid grid. 
Results were obtained for a cascade of Joukowski airfoils and are 
quite satisfactory, except near the stagnation line and the special, 
five-sided cell of the C-H hybrid grid. The unrealistic behaviour of 
the pressure near the five-sided cell was tracked down to the use of the 
SCL form of the equations. Use of the CRCL form along a single 
coordinate line containing a downstream face of the five-sided cell 
alleviated this difficulty. 

In concluding, two important observations are made regarding the 
analysis and solution procedure developed in the present study. The 
implementation of the boundary conditions in an implicit manner, that is 
compatible with the finite difference equations employed in the interior 
of the computational domain, made it possible for the present procedure 
to provide satisfactory solutions, without requiring any externally 
added artificial viscosity. Also, the procedure is capable of producing 
satisfactory solutions for compressible flows, with no modifications 
being needed for analyzing nearly incompressible flow as well. This is 
generally not true for most other available densi ty— based formulations. 

As mentioned earlier, the solution procedure was sensitive to the 
time step At. In future work, an alternating-direction implicit (ADI) 
method, rather than the present ADE type of procedure, could be employed 
to solve the semi-elliptic equations. Also, the use of a strongly- 
implicit (SI) procedure, which is known to be less sensitive to the 
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problem parameters, should be explored, in conjunction with the IPNS 
model. The present analysis can also be extended to consider unsteady 
flow problems; important practical applications, such as rotor-stator 
interactions, could then be considered. Also, turbulent flows can be 
analyzed by the inclusion of a turbulence model such as the Baldwin- 
Lomax model [54]. Turbulent flow solutions could then be compared with 
experimental results. 

The usage of the hybrid C-H grid with the semi-elliptic model has 
exposed a few other areas for further research. Firstly, for low 
subsonic flows, use of the CRCL form is quite satisfactory. It does not 
place any requirements on the representation of the metric terms, 
especially near special points and cells, where the metrics could be 
discontinuous along a particular coordinate direction. Secondly, 
decomposition of the conservation form of the governing equations of 
motion along the transformed-coordinate directions, using the 
contravar iant components of velocity instead of the Cartesian 
components, may be more appropriate. This form of the equations is more 
natural for analyzing flows in general geometries, as they are the 
fundamental conservation equations written directly in the body-fitted 
coordinate system, rather than being transformed from the Cartesian 
coordinates to this system. This interface in the hybrid grid became 
exposed as a boundary only because of the semi-elliptic formulation 
employed. It would remain as an interior computational line if the 
complete NS equations were employed. Hence, an NS analysis on hybrid 
grids is highly desired. The use of complete NS equations would also 
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resolve some of the difficulties arising due to the change in the 
orientation of the coordinate directions and the use of a coordinate- 
related approximation in the vicinity of the five-sided cells. 
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APPENDIX A 


INVERSION OF COEFFICIENT MATRIX OF BLOCK-TRIDIAGONAL 
SYSTEM WITH PERIODIC BOUNDARY CONDITIONS 

The system of equations resulting from the discretization of the 
governing equations along a periodic boundary, where the peridicity 
condition is imposed in an implicit manner, is of the following form. 


[p] AQ = D (A.1 a) 

In expanded form, Eq. (A.1 a) is written as follows. 


B 1 C 1 C ... 0 0 A 


'AQi 


1 

cT 

1 

a 2 b 2 Z p ... 0 0 0 

I 


AQ 2 


°2 

• 

j 

• « * • • 

0 0 - ... a n-1 b nh c n _ 1 

c n 0 : ... o a n b n 


4 «n-i 

“n 


* j 

°N-1 

d n 







(A.1 b) 


Here, the subscript N corresponds to the index (JMAX-1) in Eq. (3.44). 


Also, the individual entries A., B., C. are (4x4) blocks, while AQ . and 

J J J J 


D . 
J 


are (4x1) blocks. 


The inversion of the coefficient matrix [p] is performed using a L- 
U decomposition, i.e., by factoring P into upper and lower triangular 
matrices U and L, respectively, such that 
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P ' L U 


(A. 2 a) 


where 


' 1,1 


L 2,1 L 2,2 


L 3,1 L 3,2 L 3,3 


and 


%i 

I U 
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1 ,2 

I U 


2,3 ' 

I U 


3,4 


%N 


1 ,N 


2 , N 


3,N 


0 


I 


(A. 2 b , c ) 


Here, I is the (4X4) identity matrix, and L. and U. are (4x4) 

J , K J , K 


blocks. 

By forming the product LU and equating it, element by element, to 
the coefficient matrix P, the expression for the elemental blocks of L 
and U can be obtained and are given as follows. 
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L 1,1 = B 1 


4,1 = *2 


L 3,i = 4,1 " 

4,1 ” C N 

U l,2 ‘ B l' 

U l,3 - U 1.« " 
1-1 

U = B. A- 
U 1,N 1 1 


L N-t,1 * ° 


_ TJ =0 

' U 1 ,N-1 


-1 


4,2 - b 2 - 4 4 4 


4,2 = A 3 


4,2 " 4,2 


- L = 0 

4-1,2 


L = - C . B. C. , 

4.2 N 1 1 

4.3 = (B 2 ' A 2 4 1 4 ) 1 4 


4,4 ” 4,5 


- U =0 

U 2,N-1 


. - (B , - A b;’ C,)"’ A 2 B,' A, , 


U 2,N 

etc. 


Hence, the Inversion of P is carried out in two steps, 
step is a forward elimination step in which the equation 

L q = D 

is solved for q. 

This is done in the following manner. 


(A. 3) 

The first 

(A. 4) 
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For m-1 


L i,i ' B . 


Vi ■ c » 


"l.K ‘ B V 4 , 


<1, - B-' D, . 

For m equal 2 through N^I f 

L * B - a n 

m m 3-1 ,m ’ 

Q m ■ L ^ (D“AD ^ 

® m # ni m t *n m — 1 1 


u 

m,m+1 

^m, m C- 

9 

u 

m,N 

- 

- l' 1 A 
m,m - 

u 

m-1 ,N ’ 

D * 
N 

°N 

L N,m-‘ 

’ 

^N,m 

- • 

' L N , m- 1 J 

m-1 ,m 


The last unknown vector ls obtained as follows. 

l n,n-i a n l n,.v-i » 

D N " °N “ L N,*M 5 N-1 * 

L N,N “ B N ' L N,N-1 (C N-1 + U N-1,N ) » 

30 that, finally, 


% * L N,N °N 

The second step is the backward elimination 
equat ion 


step, during which the 
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(A. 5) 


U AQ = q 

. .3 From the structure of upper 

1. solved to yield the solution vector M. 

■ , 4.1 nr hiock matrices along tne 

triangular »atrix U, which has identity bloc* ta 

• .inns the superdiagonal and the 

j nnn _ zero block matrices along tne a v 

It is seen that the renewing recursive relation should be 
last column, it is seen 

employed to obtain AQ. 


. „ a q - a for, m = N-1, N-2. 
AQ m + U m,m+1 A Vl ®.H N ^ 


» 1 » 


with AQ^ = 
being used to 


initiate 


the 


determination 


of the solution vector 


AQ. 
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appendix b 


Bj PRESEXTATIQ,, of THE METRTr COEFFTripwro 

10 °° em ° lent3 - ari31ng due to the transformation of the 
8 ~ g ~ - — ««) coordinates to the 

— ationsi coordinates, should he evaluated in a manner 

~t uith the disoretitatioh of the governing equations 

The metric coefficients for a general m a „ • 

(, „ g al " lappln 8 are given in Eq. 

' • ' 3 / and are ^ 


= J y n 

x n 


Zy - - J X 

y n 


n x = ' J and n = j x . 

5 (B.i a-d; 

The metric derivatives y ariH 

n 3nd "n 3re ev ^ a ted using central 
differences at a general into • 

x) . . “ oomputational point Henoe , 

n l i + 1 >j x i + i,j + l x i + 1 1 j _ 1 )/ 2 An 

and y I . = ( v 

0 1+1,J i+1 ,j + 1 ~ y i+i f j_i )/2An . (b ^ 

“ 3 b0UMary P ° Int ' the - are evaiuated using a f lrs t o d^ 

accurate representation. Accordingly, t-at-order 

At j = JMIN, 

n| i+ i»j ^ x i + i,j+i " x i+ 1 > j)/An , 

^ y i+i,j+i -y i+] ( j) / An , 

and, at j - jmax, 

- Vi.j-i >/4n , 




(B.3 a-d) 
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In order to determine the appropriate representation of x 5 and y f . 

a simple test of uniform-flow computation is employed. This test 
involves assuming a uniform flow at a streamwise location i and applying 
the numerical integration algorithm, along with uniform-flow boundary 
conditions, to obtain the solution at station (lei). If the algorithm 
is formulated consistently, a uniform flow should result at station 

to this end, it is sufficient to consider the inviscid part of 

the governing equations, i*e., 


* * 

M + 3£ , o 

35 

since ^-derivatives appear only in the streamwise convective and 
pressure-gradient terms. The quasi-linearized form of the above 
equation can be obtained from Eqn. (3-27) and is 

A 1 A X Q + A5 (B 1 A X Q) = S 


(B.4) 


where 


i* 1 5 1 -«|s <5n 


(B . 5 a,b) 


Here, 


E * 1 - (y^ E * 1 * F * 1 ' 

«* 1 9 ‘ - & 1 *’ E * 1 * ( V 1 f * 1 


g n q 1 . i ) 1 * 1 E * 1 * (?) U1 E n 


and 


(B.6 ) 
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Following the discretization procedure, as discussed in Chapter 3> 
viscid terms, the discretized form of the right-hand side of 


(B.5 b) 

at location 

(i+1 

.j) 

can be written as 



s ‘ 

n 

-*i 

E. 

J 

* 

] Fj 1 ! 


+ 1 
4 

- 

& 

|i+ ; 

j + i 

] * kH, - 


i -*i i 

F . } 

J + 1 1 



& 

i + 1 

J-1 ' 

1 * I&j-, - 

fe' J 


+ _££ 

2 A 

n Sj* 

i 

-1 

£ 

i) 1 *' g* 1 
J j"1 j-1 





-#i 

F . 1 
J“1 

- (; 

i) i + 1 p*i i 
J Vi j-W 


(B 


in order to recover a unifor ra -fi ow solution, it will be sufficient 

t0 3h ° W th3t the S0UrCS tS ™ S is starting f rom uniform-flow 

conditions at station i, and using uniform-flow boundary conditions. To 

show this, it is recognized that, for the case wi t-h 

, i or me case with constant density, 

velocity and temperature, 

P - Pi 

u = Uj 

v = 0 
and T = T 

* i * 

Where P,, u, a„ d T, are constants. Hence, the flux veotora are constant 
With respect to the n direction, i.e., 

g*i = . =*i 

j 1 j E j + 1 1 for j = JMIN+1, ... , JMAX-1 
• for j = JMIN, .... j MAX . 


-*i 

and F = o 
J 


(B.9) 
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Kith these conditions, the source term S can be written in the 
following mannner. 

S. f I [(£)* - £)]*'] 


ff— l i + 1 - f r ' x 'l i + 1 ll 

2 An Uj i J-! I— Jj.T ]} Ej 1 


(B.10) 


Substituting for ( 5/J ) and ( n/j) fr0 „ Eq . ^ enploying ^ 

discretization given by Eq. (B.2), Eq. (B.10) can be written as 
s ' 1 FaH t - <?!»,, j„ - > r i» ll j-, > * - y^.,)] 

* rS " t * ' y i*l,j ) * f y l>J e 2 - Xj j) 


_ f 

\ 


(;/ i*i,j " , i*i,j-2 ) * ^i.j - y ltJ . 2 )] 
* rh [(y s : :H,j., - ij 1 . 


(B . 1 1 ) 

For S to be equal to zero, it is clear, from the above equation, 
that y ? is to be evaluated, as an average of its neighbours, as follows. 

y s|i + 1 ,j + 1 - TTi - n.juB * - yiJ ] 

+ _J r 

2 AC y i + 1 ,j + 1 y i, j + 1 ] (B . 1 2 ) 

Similarly, by assuming v = constant and ,, n 

5 constant and u = 0, we can arrive at the 

following representation for the term x 
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(B.13) 


i+I.J " X i» 


X c|i + 1,j + 1 4 AC ^ X i+1,j + 2 X U + 2 

+ 2 AC ^ X i + 1,j + 1 X i,j + , >^ 

The representation of x^ and at n = constant boundaries, for 
example, at j - JMAX can be obtained by forming the discretized form of 
Eq. (B.5 b) at j = JMAX-1 . 

. . r i + i i ( — JMAX-1)], the source term is written as 
At the location [i+1,J l- ' ‘ J * 

follows. 


- I l - <y i 


ns 1 ’ yi *' 

•rs' - - y x.J> 

* rh ,j*t ‘ tu ' y i*i .j-’ y i.j-'' 

“ Fat 1 y i»i ,j “ y i.J * y ‘-j“ 2 ^ 

in order to produce a zero source term S, from Eq. (B.1D. y ? term 
should be represented in the following manner: 


* (y i.j*1 " 'i.J-’’ 1 


(B.1 4) 


E, i+1 , j = J MAX 


_ r v - y. .) / AC 
' y i+1 ,j i.J 


(B.15) 


Similarly, the discretization for the term is performed as follows. 


1 = fx - x. .) / AC 

‘d i+1 , j = JMAX i + 1 , J X *J 


(B . 1 6 ) 
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APPENDIX C 

DISCRETIZATION OF DERIVATIVES AT CORNERS OF FIVE-SIDED CELLS 
In a hybrid grid, ’special-points’ commonly occur at the congruence 
of more than two different types of sub-regions. As stated in Chapter 
5, these ’special-points’, B£ and E £ in Fig. 31, are presently 

considered to be the centers of non-standard 5-sided cells. In the 
discretization of the various derivative terms in Eq. (5.1), special 
consideration is required in arriving at the appropriate neighbours for 
the vertices of these 5-sided cells, i.e., in defining the computational 
molecule at the vertices of these cells. The vertices of these non- 


standard cells are B’’ 


B 2 


B" and E’’ , 
5 1 


E 2 


E 5 


as shown 


in Fig. 31. The eight immediate neighbours for each of these points are 
shown as the encircled points in the sketches ‘below. 


At vertex B’^ , 


At vertex B£ , 



The computational moleclues shown were arrived at by examining the 
multi-block structured computational region and seeking the natural 
neighbours of a given computational point. The neighbours at the 
corners of the cell around E£ are assigned in a similar manner. At a 

given vertex, the computational molecule employed during one step of the 
two step solution procedure, for the grid-generation equations as well 
as the flow equations, is the same as that for the second step of the 
procedure. This is essential for consistency of the discretized 
equations used during the two steps. 
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FIG. 2. TYPES OF 
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u,v,T 

MMmmm 

y x 

• c,x,y,e 

• c,x,y,e 


mzmmmr 


• c,x,y,e 
c,x,y ,e 


5 

4 

3 

2 

1 


NO. OF UNKNOWNS 


24 


c x 

u,v,T (lagged) 


NO. OF EQUATIONS 

5 CONTINUITY 
4 x-MOMENTUM 

(4+1) y-MOMENTUM 
4 ENERGY 

6 BOUNDARY 
CONDITIONS 


TOTAL 24 


c. WALL-WAKE BOUNDARY CONDITIONS 


m/m 


p,u, v,T 

. 5 

• c,x,y,e 4 
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mm 
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c,x, y,e 
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NO. OF UNKNOWNS 


24 


NO. OF EQUATIONS 


5 

4 

4 

4 

7 


CONTINUITY 

X-MOMENTUM 

y-MOMENTUM 

ENERGY 

BOUNDARY 

CONDITIONS 


TOTAL 24 


d. WAKE-WAKE BOUNDARY CONDITIONS (REGION-PERIODIC GRID) 


FIG. 2 (CONCLUDED). TYPES OF BOUNDARY CONDITIONS FOR CASCADE FLOWS. 
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Inflow 

Boundary 


Outflow 

Boundary 


Fully 

Developed 

Flow 

Profile 


Uniform 

Velocity 

Profile 


FIG. 1». 



FIG. 3. TWO STEPS OF THE SOLUTION PROCEDURE. 
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b. FLATf*- PLATE CASCADE 


CONSTRICTED CHANNEL AND VARIOUS CASCADE CONFIGURATIONS, 
AND INFLOW-OUTFLOW BOUNDARY CONDITIONS. 
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Uniform 

Velocity 

Profile 
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Pressure 


c. CASCADE OF EXPONENTIAL AIRFOILS 
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Velocity 

Profile 
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Constant 

Static 

Pressure 



Uniform 

Velocity 

Profile 


d. CASCADE OF PARABOLIC- ARC AIRFOILS 



Constant 

Static 

Pressure 



e. CASCADE OF JOUKOWSKI AIRFOILS WITH MODIFIED LEADING EDGES 


FIG. 4 (CONCLUDED). CONSTRICTED CHANNEL AND VARIOUS CASCADE 

CONFIGURATIONS, AND INFLOW-OUTFLOW BOUNDARY CONDITIONS. 
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0.0 



si . t/h ■ 0*1 

FIG. 6. WALL-PRESSURE AND WALL-SHEAR PARAMETER DISTRIBUTIONS 
FOR EXPONENTIAL CHANNEL CONFIGURATION, Re - 1500. 


116 






?q 29 X 30 31 

,- a x 75 100 27 28 C 

enlarged horizontal scale 




50 X 75 


100 27 26 29 x 30 31 0 

ENLARGED horizontal SCALE 


c. t/h “ 0.2 


6 (CONCLUDED). WALL-PRESSURE AND WALL-SHEAR PARAMETER DISTRIBUTIONS 
FOR EXPONENTIAL CHANNEL CONFIGURATION, Re - 1500. 
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ENLARGED HORIZONTAL 
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ENLARGED HORIZONTAL 


b. Re - 3100 


VELOCm Am 

°F FINITE put P, * 1BUT10 " S roB * CASCADE 
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enlarged horizontal 



75 9 



ENLARGED HORIZONTAL £ 


c. Re - 6200 





25 X 


ENLARGED horizontal scale 




to X 


ENLARGED HORIZONTAL scale 


d. Be • "000 

will-pressure, hake-centerline velocity 
FI0 ' 8 Tl-shear paraheter distributions for a cascade 
of FINITE FLAT PLATES. 


ENLARGED HORIZONTAL SCALE 




ENLARGED HORIZONTAL SCALE 


e. Re - 16000 


FIG. 8 (CONCLUDED). WALL- PRESSURE , WAKE-CENTERLINE VELOCITY AND 
WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF FINITE FLAT PLATES. 






5x10 



PARAMETER DISTRIBUTIONS FOR VARIOUS Re FOR A 
CASCADE OF FINITE FLAT PLATES. 
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enlarged horizontal scale 



75 8 




(iv) 

i 

- — 1 1 1 

J A 1 


" 10 X n , 

enlarged horizontal scale 


a. t/c - 0.05 


FIG. ,0. WALL-PRESSURE, WAKE-CENTERLINE VELOCITT AND 

WALL-SHEAR PARAMETER DISTRIBUTIONS TOR A CASCADE 
OF EXPONENTIAL AIRFOILS, Re . ,BOO. 






enlarged horizontal scale 




enlarged horizontal scale 


b. t/c - 0.075 


,0 (COR'D). HALL-PRESSURE. WAKE-CENTERLINE VERITY AND 
WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF EXPONENTIAL AIRFOILS, Re - 1500. 
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FIG. 10 (CONCLUDED). WALL-PRESSt 
WALL-SHEAR PARAMETER DI* 
OF EXPONENTIAL AIRFOILS, 
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PARAMETER DISTRIBUTIONS FOR VARIOUS t/c FOR A 
CASCADE OF EXPONENTIAL AIRFOILS. 




ENLARGED HORIZONTAL SCALE 
a. Re ■ 3100 

FIG. 12. WALL^PRESSURE, WAKE-CENTERLINE VELOCITY AND 

WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF EXPONENTIAL AIRFOILS, t/c - 0.05. 








ENLARGED HORIZONTAL SCALE 


b. Re ■ 6200 

FIG. 12 (CONT’D). WALL- PRESSURE, WAKE-CENTERLINE VELOCITY AND 
WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF EXPONENTIAL AIRFOILS, t/c - 0.05. 








ENLARGED HORIZONTAL SCALE 
c. Re - 11000 

FIG. 12 (CONT'D). WALL- PRESSURE , WAKE-CENTERLINE VELOCITY AND 
WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF EXPONENTIAL AIRFOILS, t/c - 0.05. 
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FIG. 


ENLARGED HORIZONTAL SCALE 
d. Re ■ 15000 


2 (CONCLUDED). WALL^PRESSURE , WAKE-CENTERLINE VELOCITY AND 
WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF EXPONENTIAL AIRFOILS, t/c - 0.05. 


3 








135 


° (186x71) 



EXPONENTIAL airfoils, t/c 
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X 


FIG. 1 


ENLARGED HORIZONTAL SCALE 


a. t/c - 0.05 


5. WALL- PRESSURE, WAKE-CENTERLINE VELOCITY AND 

WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF PARABOLIC-ARC AIRFOILS, Re - 1500. 







ENLARGED HORIZONTAL SCALE 


b. t/c - 0.075 

FIG. 15 (CONT'D). WALL* PRESSURE, WAKE-CENTERLINE VELOCITY AND 
WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF PARABOLIC- ARC AIRFOILS, Re - 1500. 
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FIG. 


ENLARGED HORIZONTAL SCALE 


c. t/c - 0.1 


5 (CONCLUDED). WALL-PRESSURE, WAKE-CENTERLINE VELOCITY AND 
WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF PARABOLIC-ARC AIRFOILS, Re - 1500. 
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^rrr™ ™ ind 

pababolic-jrc airfoils, t/o . 0.05. 
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, u»i.ly- PRESSURE, *AKE*CE*TERLINE VELOCm RED 
F10, 17 rJs^R PARAMETER DISTRIBUTES R. A CASCADE 

0F p »rabolic-arc airfoils, t/e - 0.05. 
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ENLARGED HORIZONTAL SCALE 




d. Re - 15000 


FIG. 17 (CONCLUDED). WALL- PRESSURE , WAKE-CENTERLINE VELOCITY AND 
WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF PARABOLIC-ARC AIRFOILS, t/c - 0.05. 
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PARAMETER DISTRIBUTIONS FOR VARIOUS Re FOR A 
CASCADE OF PARABOLIC*- ARC AIRFOILS. 
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enlarged horizontal scale 
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enlarged horizontal scale 


a. M - 0.032 


FIG. 19. WALLfPRESSURE , WAKEfcCENTERLINE VELOCITY AND 

WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF PARABOLIC* 4 ARC AIRFOILS, t/c - 0.05, Re 


ENLARGED HORIZONTAL SCALE 




enlarged horizontal scale 


b. M - 0.128 


FIG. 19 (CONT'D). WALL- PRESSURE, WAKE*CENTERLINE VELOCITY AND 
WALL*- SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF PARABOLIC*’ ARC AIRFOILS, t/c - 0.05, Re - 15000. 




ENLARGED HORIZONTAL SCALE 


c. M - 0.M9 
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FIG. 19 (CONCLUDED). WALLr PRESSURE, WAKE-CENTERLINE VELOCITY AND 
WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF PARABOLIC-ARC AIRFOILS, t/c - 0.05, Re - 15000. 









800 * 



CASCADE OF PARABOLIC'iARC AIRFOILS 









NO. VALUE 

1 0.29226E 01 

2 0.29313E 01 

3 0.29900E 01 
9 0.29987E 01 

5 0.29579E 01 

6 0.29662E 01 

7 0.29799E 01 

8 0.29836E 01 

9 0.29923E 01 

10 0.3001 IE 01 

11 0.30098E 01 

12 0.30185E 01 

13 0.30272E 01 
19 0.30359E 01 

15 0.30997E 01 

16 0.30539E 01 

17 0.30621E 01 

18 0.30708E 01 

19 0 . 30796E 01 

20 0.30883E 01 

21 0.30970E 01 


FIG. 21 . STATIC PRESSURE CONTOURS FOR A CASCADE OF 

PARABOLIC-- ARC AIRFOILS, Re - 15000, M - 0.49. 

QO 
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NO. VALUE 

1 0.29226E 0 

2 0.29313E 0 

3 0.29900E 0 
9 0.29987E 0 

5 0.29579E 0 

6 0.29662E 0 

7 0 • 297*f9E 0 

8 0.29836E 0 

9 0.29923E 0 

10 0.3001 IE 0 

11 0.30098E 0 

12 0.30185E 0 

13 0 . 30272E 0 
19 0.30359E 0 

15 0 • 30997E 0 

16 0.30539E 0 

17 0.30621E 0 

18 0.30708E 0 

19 0.30796E 0 

20 0.30883E 0 

21 0.30970E 0 


FIG. 21b. ENLARGED VIEW NEAR LEADING EDGE REGION. 
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FIG. 23. WALL-PRESSURE, WAKEr-C ENTERLINE VELOCITY AND 

WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF JOUKOWSKI AIRFOILS (MODIFIED LEADING EDGE). 






FIG. 23 (CONCLUDED). WALL- PRESSURE, WAKE-CENTERLINE VELOCITY AND 
WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF JOUKOWSKI AIRFOILS (MODIFIED LEADING EDGE). 
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FIG. 2M. CONVERGENCE HISTORY FOR A CASCADE OF PARABOLIC-ARC 
AIRFOILS, t/c - 0.05, Re - 6300. 




FLAT- PLATE CASCADE WITH Ax - CONSTANT. 








- 1 



0 1 00 200 300 

No. of Iterations 


FIG. 27. EFFECT OF At ON CONVERGENCE FOR A CASCADE 
OF EXPONENTIAL AIRFOILS. 
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28. WALL-PRESSURE, WAKE-CENTERLINE VELOCITY AND 

WALL-SHEAR PARAMETER DISTRIBUTIONS FOR A CASCADE 
OF EXPONENTIAL AIRFOILS, t/c - 0.05, Re - 25000. 





a. CONFORMAL C-GRID (AFTER REF. 41) 



b. * LINE-PERIODIC * H-GRID 


c. ’REGION-PERIODIC 1 H-GRID 


FIG. 29. VARIOUS TYPES OF GRIDS FOR CASCADES. 
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a. OVERSET GRID 


patch 

BOUNDARIES 



PATCH 

BOUNDARY 


ZONE 3 


PATCH 

boundary 


b. PATCHED-DIS JOINT GRID 



c- PATCHED’- JOINT GRID 


FIG. 30. 


TOES OF MULTI-BLOCK STRUCTURED GRIDS 
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a. CELL WITH NON-STANDARD NUMBER OF FACES 



b. POINT WITH NON-STANDARD NUMBER OF NEIGHBOURS 


FIG. 33 . TYPES OF SPECIAL CELLS AND POINTS. 
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a. STAGGERED CASCADE (GOSTELOWS CASCADE) 



U..STAGCEEED CASCADE OE JOUKOWSKI AIHFOILS 


FIG. 37. 


typical C-H HYBRID GRIDS. 


m 
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b. 


a. TWO-STEPS OF THE SOLUTION 
PROCEDURE 


b. ENLARGED VIEW OF LE REGIONS 



HATCHING DIRECTIONS FOR SEMI-ELLIPTIC FLOW- 
SOLUTION PROCEDURE ON C„H HTBRID oniD . 


172 


0.8 


0.20 



a. SURFACE 


PRESSURE AND WAKE-CENTERLINE VEIDCITY DISTRIBUTIONS 



FIG. 39. RESULTS FOR A CASCADE OF JOUKOWSKI AIRFOILS 
USING C-H HYBRID GRID, Re ■ 150. 
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PRESSURE 


“S FOR A C 
GRID, Re 


i 



d. STATIC PRESSURE CONTOURS SUPERIMPOSED ON THE HYBRID GRID. 


FIG. 39 (CONCLUDED). RESULTS FOR A CASCADE OF JOUKOWSKI AIRFOILS 
USING C-H HYBRID GRID, Re - 150. 
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a. SURFACE PRESSURE AND WAKE-CENTERLINE VELOCITY DISTRIBUTIONS 



b. WALL- SHEAR PARAMETER DISTRIBUTION 


FIG. 40. RESULTS FOR A CASCADE OF JOUKOWSKI AIRFOILS 
USING C-H HYBRID GRID, Re - 600. 
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d. STATIC PRESSURE CONTOURS SUPERIMPOSED ON THE HYBRID GRID 


FIG. 42 (CONCLUDED). RESULTS FOR A CASCADE OF JOUKOWSKI AIRFOILS 
USING C*H HYBRID GRID AND COMBINED SCL-CRCL 
FORM OF EQUATIONS, Re - 300. 
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decomposition of the physical domain leads to^multi ^lock so i ut ion procedure *hi ch uni i ke 

physical -problem boundaries. This P*™ nrocedure Satisfactory results are obtained for flows 
most found in literature, is not a pateh^g p 0 ■ „ the i PNS formulation on the C-H gnd 

through cascades of Joukowski airfoils. ep require special treatment. However, with 

exposes two small portions of the ^d interfaces a Recommended, so as also to avoid mcon- 
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